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ABSTRACT 


The considerations in selecting- the sampling rate for a digital 
control of aircrafts are identified and evaluated. The design method 
used in this analysis is the optimal discrete synthesis. Due to dis- 
cretization of the continuous plant this method of design does not 
introduce an artificial limitation on the sampling rate. The principal 
example used is the short period mode of a hypothetical high perform- 
ance aircraft. The assumed model includes a bending mode and wind gusts. 

Four major factors which influence the selection of the sampling 
rate are identified: (a) the time response to control inputs; (b) the 

response to an external disturbance; (c) the sensitivity to variation 
of parameters; ( d) the roughness of the response to control inputs. 

Each of these factors and its relation to the sampling rate was inves- 
tigated. It was found that the limiting factors in the selection of 
the sampling rate for the example are the time response to a control 
input, and the response to an external disturbance. The sensitivity to 
variation of parameters is larger for lower sampling rates. However, the 
sensitivity can be reduced by modifying the design of the optimal linear 
compensator. Different roughness functions which measure the roughness 
to control inputs are suggested and demonstrated on the example. 

The optimal discrete synthesis computer program, which is based on 
eigenvector decomposition of the state-costate Hamiltonian matrix, is 
a highly efficient program. This program calculates the optimal discrete 
regulator, the steady state Kalman filter, and the mean response to 
external disturbances. 
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I . INTRODUCTION 


A- PROBLEM STATEMENT 

During the last ten years, the aerospace industry has shown a ten- 
dency to replace the analog components of closed loop systems with digital 
computers # The advantages of using a digital computer instead of a 
specially built analog system are numerous. Among them are greater 
accuracy, the ability to change parameters during the control operation, 
and flexibility of the control logic. On the other hand, the principal 
disadvantage of a digital computer in a closed loop control system is 
its discrete mode of operation. The computer processes numbers generated 
in real time by sampling continuous signals. The computer outputs, which 
are sequences of numbers, have to be reconstructed into analog signals 
(commands to actuators). Therefore, in the process of designing a digi- 
tal autopilot, careful selection of the rate of sampling and the process- 
ing of commands are important. 

Selecting an appropriate sampling rate for an aircraft digital 

controller necessitates a compromise. Cost and accuracy are factors 

which argue for lowering the rate, go . A low co directly reduces 

s s 

the cost of A/D and D/A equipment. Using less central processing unit 
percentage time can either free the system for other functions, or 
result in reduced central processor costs. The increased accuracy ob- 
tained by slower sampling is well documented [BO-1] and can be trans- 
lated into additional cost savings by reducing the word size. Economic- 
ally speaking, the best engineering choice is the slowest possible samp- 
ling rate meeting all performance specifications. 

Factors which may constitute an incentive to increase the sampling rate 
are, e.g. , (1) closed-loop bandwidth, or time response requirements; 

(2) sensitivity to parameter variations; (3) effect of random disturb- 
ances; and (4) roughness of control. 
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Our objective will be to identify and analyze the important factors 
which influence the selection of a sampling rate for a closed-loop air- 
craft discrete control system. After discussing the factors and their 
properties, we will propose methods which will help the designer choose 
a sampling rate which will not violate given criteria of performance 
or given properties of the system. 

As an example, we chose a hypothetical, Mach 3 aircraft [ BO-1, SU-l], 
flying in a highly turbulent atmosphere. This choice, which we shall 
term F-H , was made because the requirements of a high-performance military 
aircraft impose a limit on the minimal sampling rate. 

B. PREVIOUS RELATED RESULTS 

Aircraft digital control systems have now been implemented [DE-1, 
MA-1], and many more have been studied [BE-1, BO-1, DS-1, LE-1, SU-1], 

The pertinent literature discusses various methods for selecting the 
sampling rate. 

Sampling rate selection is sometimes based on a specific multiple 
of the highest important bending mode. An appropriate value for this 
multiple was reported by Lee [LE-l] to be four. Others [JO-1, SI-1, ST-1] 
have also selected sampling rate to be about four times the highest 
important bending mode; however, it is perhaps more typical to select 
sample rates at approximately six to ten times the highest bending mode 
[see, for example, Refs. BO-1, DE-1, ED-1, OS-1, SU-1], Recently, 

Berman [BE— 1] introduced the concept that the proper sampling rate is 
independent of the bending modes and should be based solely on disturb- 
ance effects. This was applied to a V/STOL example and yielded a samp- 
ling rate which was slower than the highest bending mode. According to 
some discussions [LE-1], such a slow sampling rate violates the sampling 
theorem [RA-1], which is said to state that the sampling rate must be 
at least twice the highest bending mode frequency. Our interpretation 
of the sampling theorem will be discussed in Chapter IV, 

A completely different approach to selection of sampling rates 
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was taken by Mrazek [MR-1], His selection process was dominated by the 
time constants of the measuring 1 instruments 1 analog prefilters. This 
consideration resulted in sampling rates 100 to 200 cps, However, 
flying experience with digital autopilots [MA-1] shows that such high 
sampling rates as suggested by Mrazek are unnecessary. 

The design technique used in the analysis often influences the samp- 
ling rate selection and, in some cases, causes the sampling rate to be 
significantly faster than required. Although many variations exist, 
we prefer to divide the design methods into two broad categories: 

(1) those whose design is done in the continuous domain (or s-plane), and 

(2) those whose design is done in the discrete domain (z or w~plane). 
Designs of the first category are attractive since they utilize the 
experience gained over many years of continuous autopilot designing. 

The authors using this method [BO-1, ED— 1 , OS— 1, SU-l] choose one or 
a combination of discrete approximation techniques (reviewed recently 
by Edwards [ED-2] and Slater [SL-1]) to transform the resulting contin- 
uous compensation into a discrete compensation. The effect of the ap- 
proximation in the design is typically checked by a precise simulation. 

Designs of the second category include the w-plane techniques 
[LE-1, ST-1], z-plane Nyquist techniques [SI-1], and the discrete state 
space techniques of Berman [BE-1], Johnson [JO-1], and the author of 
this work. 

The approximations inherent in category 1 (s-plane) methods introduce 
an additional constraint which may be important in sampling rate selection. 
It is interesting to note that all authors reporting the use of a cate- 
gory 1 design method selected a sample rate which was a higher multiple 
of the bending modes of interest than those authors using a category 
2 method. 

An attempt was made to relate the sampling rate to overall per- 
formance for the case of optimally controlled systems. Lewis and 
Athans [LEW-1], and Astrom [AS-1] developed a method which investigated 
the change in the quadratic index of the continuous plant for different 
sampling intervals. Their results, which are essentially experimental 
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(on computer), show that for larger sampling intervals, the general 
tendency of the quadratic index is to increase* 


£. NEW RESULTS 

Four major factors which influence the selection of the sampling 
rate for the discrete control of a continuous system were identified 
and analyzed. These factors are: 

1. the time response, 

2. the response to an external noise, 

3* the sensitivity to variations of parameters, 

4. the roughness of control. 

By using an optimal discrete design, we eliminated the time lags 
introduced into the control loop by the discretization of a continuous 
design. Hence, the sampling rate is no longer limited by discretization 
approximations and we can concentrate on the analysis of the four factors 
listed above. 

1. The time response to a step input is directly related to the 
length of the sampling interval. For a given criterion in the 
time domain, the time response deteriorates for longer sampling 
intervals. From simulation of the F-H short period mode, we 
found that the discretized low pass filter and the discretiza- 
tion of the pilot's analog input introduce a considerable delay. 
This delay is as long as two sampling intervals and limits the 
choice of the sampling rate. For the F-H short period example, 
the lower limit of the sampling rate, imposed by the time re- 
sponse, was found to be in the vicinity of 10 cps. 

2. The response to external disturbances is a function of closed 
loop dynamics, the disturbance correlation time, and the samp- 
ling interval. It was found that for any available (discrete) 
control, there is a well-defined limit to noise alleviation. 

This limit is a function of the sampling rate. For higher 
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sampling rates, the noise alleviation by the closed loop system 
is greater. 

By analyzing the external disturbance as a Gauss-Markov 
process, and by using an optimal discrete control, greater 
noise alleviation for the same sampling rates can be achieved. 

3. Sensitivity to variation of parameters was reduced by proper 
modeling of the discrete compensator. The optimal discrete 
compensator assumes a perfect knowledge of the system and the 
noise. However, if some of the system’s parameters vary from 
their nominal value, the filtering and controling action of the 
optimal compensator is distorted. This distortion was shown to 
be worse at lower sampling rates. 

In the F-H short period example, a 10 percent uncertainty in 
the bending frequency causes instability of the closed loop 
system. By remodeling the discrete observer, essentially by 
increasing the damping of the observer’s error poles corres- 
ponding to the bending mode, the closed loop system was sta- 
bilized. The same could be done for larger sampling intervals. 

Further restriction in sensitivity to variations of the bending 
frequency can be achieved by a damping augmentation of the 
bending mode. 

4. Roughness of control is caused by the quantization of the input 
signal by the zero order reconstruction hold. Intuitively, it 
is obvious that for very high sampling rates, this roughness 

is negligible. But for lower sampling rates, the phenomenon 
cannot be ignored. This is more than a theoretical speculation. 

It was also reported by Mathew [MA-1], who detected an undesir-* 
able jittery action on the digitally-controlled Saab actuators. 

To put the concept on an analytical basis, different roughness 
functions (RF) were defined* Basically, the RF is defined as a 
sum (for an impulse response) ; or as a mean value (for rms response) 
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of the squares of the control discontinuities. Algorithms for 
calculating the RF are given and the RF concept is demonstrated 
for the F-H short period mode configuration. An interesting 
result came out: if the quadratic cost function is kept con- 

stant, the roughness of control may decrease for larger sampling 
intervals. This phenomenon is fully explained in Chapter VII, 

D. THESIS OUTLINE 

In Chapter II, various discrete control techniques will be surveyed. 

In Chapter III, the fundamentals of a continuous system’s discretiza- 
tion and of the optimal discrete theory will be outlined. The author’s 
approach to an optimal discrete synthesis, by eigenvector decomposition, 
will be given. Various numerical algorithms, useful in discrete analysis, 
will be described. 

In Chapter IV, the limit on the sampling rate, imposed by a time 
response, will be investigated. An example, which illustrates the basic 
concepts and is used in other chapters, will be described. 

In Chapter V, we will investigate the relation between the sampling 
rate and the response to an external disturbance, A theoretical proof 
of the limit of noise reduction as a function of the sampling rate will 
be given. 

In Chapter VI, the sensitivity to variations of parameters will be 
described. We will investigate the behavior of a closed loop system 
which includes unwanted frequencies, in our case, the bending mode. We 
will show how sensitive the system is to an imperfect knowledge of un- 
wanted frequencies and what can be done to reduce this sensitivity. 

An interesting relation between sensitivity and the stability of 
the compensator will conclude this chapter. 

In Chapter VII, a new criterion (the roughness function), will be 
explained and methods for calculating the function will be given. This 
new concept will be demonstrated in the design of the control loop of 
the F-H short period mode example. 
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In Appendix A, an instruction manual and a program DISC are given. 
This program, for the synthesis of an optimal regulator and an optimal 
steady state observer for discrete linear systems, has various options. 
The discretization procedure and the calculation of the response to an 
external disturbance are included. An illustrative example is given at 
the end of Appendix A. 

In Appendix B, the simulation algorithm of the principal example's 
behavior will be explained. This simulation, based on a specially built 
computer program, simulates the behavior of a discretely controlled con- 
tinuous system. An imperfect knowledge of various parameters of the 
simulated system is included as an option in the program. Most of our 
results regarding the time response of our principal example are de- 
duced from this simulation. 


E. SUMMARY OF CONTRIBUTIONS 

1. The first investigation of sampling rates, considering all 
effects and identifying what is important. 

2. The factors which influence the selection of the sampling rate 
were evaluated for an F-H example. It was shown that a proper time 
response and gust alleviation were the two major factors determining the 
sampling rate. 

3. Demonstration of a technique that eliminates bending modes and 
uncertainty of their exact frequencies as an important constraint on many 
aircraft autopilot designs. All cases are easier than that of an F-H 
flying at Mach 1.2 and, at zero altitude. 

4. A relation between the stability of the compensator and the 
sensitivity was described. 

5. It was found that there is a definite limit on the alleviation 
of an external noise by a discrete controller. This limit is a function 
of the sampling rate. 

6. A new criterion, the roughness function, was defined. Methods 
for calculating the function by using existing algorithms are given and 



demonstrated for the F-H example 


7. A new algorithm is given for calculating the Liapunov equation 
using eigenvector decomposition, 

8, A computer program was developed for the synthesis of an optimal 
discrete regulator and an optimal discrete steady state observer based 

on eigenvector decomposition. 



II . SURVE Y OF TECHNIQUES FOR DESIGNING DISCRETE 
CONTROLLERS FOR CONTINUOUS SYSTEMS 


The method of design, divided in the previous chapter into two 
major categories, can be further divided into the following cate- 
gories : 

A. A classical continuous design and a discretization, 

B. A continuous design and an optimal discretization, 

C. A classical discrete design, 

D. Discretization and optimal discrete design, 

_A, CLASSICAL CONTINUOUS DESIGN AND DISCRETIZATION 

Design on the s— plane and discretization of the compensation network 
is a widely used method. Borow et al [BO-l], Edwards [ED-1], Osder [OS— 1], 
and Sutton et al [SU-1] use this method of design. Their fundamental 
approach is to duplicate the analog filters by digital filters. 

Edwards [ED-2] and Slater [SL-1] investigated various methods which are 
used to design digital filters having properties similar to the corresponding 
analog filters. The transformations used in converting the filters from the 
s-plane to the z-plane are: the standard z-transform, the bilinear 

transformation, z-forms, and the matched z-transform. These transforma- 
tions yield digital filters of the same order as the original analog 
filter. Edward’s conclusion is that the z-matched transform has proved 
to be a good technique for generating all the standard filter forms. 

McGough [McG-l] analyzed this design method and suggested using the fre- 
quency response characteristics of a zero order reconstruction hold to 
filter out the unwanted bending frequencies, which usually consist of 
the highest frequency component. We will show in Chapter VI that in 
practice, sampling at the bending frequency may cause instability. 
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B. CONTINUOUS DESIGN AND AN OPTIMAL DISCRETIZATION 


Designing a discrete controller in the continuous domain is attrac- 
tive to any designer because the required characteristics and the closed 
loop properties are usually given in the time domain (or s-plane). 

Most designers have a better understanding of the physical quantities 
expressed in the s~plane. Hence, several attempts were made to trans- 
form the continuous design to a discrete domain in some optimal way* 

Melzer and Kuo [ME— 1] use Taylors approximation for the solution of 
an optimal regulator as a function of the sampling interval* Once a 
sampling interval is chosen and an optimal continuous design is done, 
the feedback matrix is obtained by Taylor 1 s series approximation. As 
the authors claim, this method can be verified only by numerical experi- 
ments. The method was later improved by Kuo and Peterson [KU-l]. 

Using Taylor 1 s expansions of the feedback matrix and the solution of the 
Riccati equation, they modified the feedback gains so that the response 
of the sampled-data model was as close to that of the original continuous 
system as possible. 

Yet the most promising approach of this kind is given by Yackel, 

Kuo, and Singh [YA-1] # Their method, based on Kuo’s previous results, is 
a complete digital redesign of continuous systems by matching states in 
multiple sampling periods. Their method is essentially based on the con- 
trollability theorem, which states that an n-order discrete system can be 
brought to an arbitrary state in no more than n-steps. The desired 
state is the solution of the continuous system during the interval n x T, 
(T is the sampling interval.) The main disadvantage of their method lies 
in the fact that the gains have to be changed for every sampling period. 
This adds to the complexity of the numerical calculation in the real 
time computer. 


£. CLASSICAL DISCRETE DESIGN 

This is an exact method based on z-transf ormation of the continuous 
plant F(s). The plant F(s) is transformed to a discrete plant F(z) 
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(via a zero order hold) on the z-plane, see Fig. II-l. F(z) includes 
the zero order hold (ZOH). 



Methods of design on the z-plane are well documented [e.g. , CA-1, 

CA— 2, RA-1 ]. However, designers who are accustomed to s-plane formula- 
tions prefer to transform the F(z) plant to the w-plane. The F(z) 
plant is transformed to a w-plane by the bilinear transformation 
z = (l+w)/(l-w). 

F(w) has many properties similar to the continuous plant, F(s) 

[e.g., R A— 1 ] . Therefore, a compensator design can be done by classical 
methods (Bode plot, root locus). 

After determining a proper compensation D(w), D(w) is transformed 
back to the z-plane, w = (z-l)/(z+l). D(z) immediately gives the re- 
cursive compensating equations. 

Analytical expressions for J[F(s)] are only known for very simple 
transfer functions. However, computer algorithms which transform F(s) 
to F(z) are available. The main disadvantage of this method is its 
inability to design properly multi-input, multi-output systems. 


_D. DISCRETIZATION AND OPTIMAL DISCRETE DESIGN 

The optimal discrete design is based on fundamental work done by 
Kalman [KA-1], He pointed out the equivalence between the optimal 
linear discrete regulator and the optimal linear discrete observer, 
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and proposed algorithms for computing them. This is an exact method since 
no discretization approximations are made. Franklin [GU-1] has shown 
that the optimal discrete compensator can be designed by first design- 
ing the optimal controller and then the optimal observer. We can 
separate this approach into a set of two distinct problems: the trans- 

formation of a continuous system into an equivalent discrete system, 
and the numerical solution of the matrix Riccati difference equation. 

The discretization of the continuous system is done by a numerical 
calculation of the state transition matrix and its integral. This subject, 
seemingly elementary, has not been exhausted. New efficient algorithms 
are discovered and rediscovered all the time. For example, Hansen jHA-1] 
rediscovered the algorithm devised nearly 25 years ago by Frame [FRA-l] 
and Fadeev. 

Numerical methods for calculation of the discrete regulator and the 
steady state filter are extensively documented [e.g., BR-1, and KW-1] . 

Those methods are based on recursive computation of the matrix Riccati 
difference equation until a steady state solution is reached. Following 
Potter [PO-l], who solved the matrix Riccati differential equation by a 
nonrecursive method, Vaughan [VA-l] found a nonrecursive solution for the 
matrix Riccati difference equation. His method involves the calculation 
of the eigenvalues and eigenvectors of the canonical state-costate equa- 
tions, Using the QR transformat ion--a highly efficient algorithm for cal- 
culating eigenvectors [FR-l], Bryson and Hall [BR-2] constructed a com- 
puter program for steady-state optimal control and filter synthesis of 
a continuous system. Independently of Vaughan 1 s results, the author of 
this work constructed a method and a computer program for a steady state 
optimal discrete control and filter synthesis by eigenvector decomposi- 
tion. 


The Vaughan solution is further explored by Howerton [HO-l], who 
shows that further simplification of the discrete algebraic Riccati 
equation can be achieved by transforming the system to a Luenberger 
canonical form. 
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III. OPTIMAL DISCRETE SYNTHESIS 


There are different methods of designing a discrete compensator for 
a continuous system. These methods can be divided into two basic cate- 
gories as explained in Chapter I. They are: (1) those in which the 

design is done in the continuous domain/s-plane ; (2) those in which the 

design is done in the discrete doraain/z or w-plane. 

The methods of the second category result in no artificial sampling 
rate constraints due to discretization procedure. Therefore, these 
methods are more suitable for our major objective in this work; i.e., 
to investigate the various factors which influence the sampling rate 
selection. The author of this work preferred the discrete state approach 
(instead of the w-plane approach) for the following reasons: (a) the 

state space approach easily handles multi-input multi-output systems; 

(b) the w-plane design method, essentially a classical method of design, 
does not have convenient computational means for the minimization of 
disturbance influences. 

It will be shown in further chapters that the optimal discrete 
approach, which is used in this work, is a valuable design tool. 

In this chapter we will describe the procedure for synthesis of an 
optimal linear regulator based on minimization of a quadratic cost func- 
tion for a linear, t ime- invar i ant discrete system. The procedure for 
synthesis of a discrete optimal filter will be described by using the 
Kalman analogy between an optimal regulator and an optimal steady state 
filter. 

The calculation of the regulator and the filter is based on eigen- 
vector decomposition of the related state-costate Hamiltonian matrix, 

A discretization procedure will be described for the case in which 
the controlled system is continuous and the control input, generated by 
a digital computer, is reconstructed by a reconstruction hold. 
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Various useful numerical algorithms will be given 


A. DISCRETIZATION PROCEDURE 
A-l Reconstruction Holds 

Most of the practical cases of designing a discrete controller are 
those in which the controlled system is an analog plant. A special hard- 
ware element, called the "reconstruction hold", is inserted between the 
computer and the controlled plant. The purpose of the reconstruction 
hold is to convert a sequence of numbers, usually equally spaced in time, 
to a continuous signal. See Fig. Ill— 1. 


u 

o 


DIGITAL 

COMPUTER 

l 

RECONSTRUC- 
TION HOLD 

u 

ANALOG 

PLANT 








5 

i \ 


y 


FIG. III-l DISCRETE CONTROL OF AN ANALOG PLANT. 

u ss control signal, Uj_ - sequence of 
numbers, T = sampling interval. 


For real time in closed loop control systems, the reconstruction hold 
generates signal u(t) (for t £ nT), based on u for i £ n. 
Basically, the reconstruction scheme is an extrapolation. The zero order 
hold (ZOH) , generates a constant signal u(t) = for iT £ t < (i+l)T. 
This reconstruction hold is widely used. Higher order holds are used 
primarily as examples in literature. The first order hold (FOH) is 
based on the last two control vectors — u. , and u, and generates a 

signal varying linearly with time: 
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u(t) 


u. + 
X 


T 


T 


U. 


i-1 


(3.1) 

iT <: t < ( i+l)T 0 £ T ^ T . 

This formulation of the first order reconstruction hold originates in 
sample data theory. The sample data approach usually considers a single 
input control channel and sequentially processed scalar quantities. 

If a state space approach and a digital computer are used, the FOH can 
be reformulated. The new FOH can be defined as 


u( t) £ + k^T 

IT ^ t ^ (i+l)T 0 £ t £ T 


(3.2) 


where and are quantities calculated simultaneously in the digital 

computer. Recall that the FOH is a hardware device, which generates 

the step and the linear rate k^ of the continuous control signal 

u(t), during the interval T. The vector k^ (a scalar for a single 

input system), can be calculated with respect to some specific criteria 

and does not necessarily have to be equal to (u.-u )/t. 

i i-l'' 


A-2 Formulation and Algorithms For Discretization Procedure 

_a. The transition matrix and its first integral . 

The discrete formulation of a linear continuous, differential 
system is essentially the solution in the time domain from sample point 
to sample point. The solution of 

x = Fx + G^u (3.3) 


is given by [KW-l] 


x 


i+1 


( i+l ) T 

®(t)x + J <&[(i+l)T 

iT 


T ]G^u(f)dx . 


(3.4) 


- 15 - 


For a time invariant system, the matrices F and are constants; 

therefore, (3.4) may be rewritten as 


r 

x = ®(T)x i + 1 ®(t)G^u(t ) dT 


(3.5) 


where the transition matrix < 5 (t) is given as 

$(t) = e FT = I + Ft + £57- + ... . 


(3.6) 


If li ( x ) is a constant during the interval T, then the second term on 
the right hand side of (3.5) is defined as r^(T)u , where T^(T) is 


T 1 (T) = J $(T)G 1 dt . 


(3.7) 


_b, Discretization of a continuous system driven by a white noise , 
The disturbed system: 


x = Fx + G w [w -> N(0, Q) ] 


(3*8) 


may be represented at sampling points [BR-1] by: 


x. = $x. + T w w N(0 Q ) 

i+l 1 2d, d . ' d 

L X 


(3,9) 


where 


r T 

Q d = \ $(OQ$ (t )dt 


o 

T 


- [ 0 (t )GdT , 

^o 


^ We will use the following notation: w -*N(0,Q), where w is a random 
disturbance, normally distributed, with a zero mean and a power spec- 
tral density matrix Q. 

* The reader's attention is directed to the fact that the symbol T is 
used for the sampling interval and also for a transposition of matrices. 
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We are interested in the covariance matrix Q d for a statistical anal- 
ysis of system behavior, whereas, and would only be needed for a 
complete stochastic simulation. Therefore, (3.9) will be replaced by 


i+1 


I = 


+ Iw , 

1 d , 


1 0 

0 1 


w 


d. ^ N(0, Q d ) 


Q = ( <Mt)G QG T 0 T (T)dT . 
CL ,) 2 


(3.10) 


c_. Discretization of a continuous cost function . 

A continuous quadratic cost function J(x, u, t): 


J = ^ (x T Ax + u T Bu)dT 



A 

A 


x . 

N-l T T 

11 

12 


i 

i“0 

A 

A 


u . 


21 

22 


i 


(3.11) 


can be transformed to a discrete version 

A i 2 i [*±1 

(3.12) 

where N = t/T; by the following procedure, (3.11) can be rewritten as 

N-l (itl)T 


J = 


= ^ f (x T Ax + u T Bu)dT . 

i=o L 


(3.13) 


iT 


The interval 0 - t was subdivided into N intervals T. The integrals 

inside the summation expression of (3.13) can be expressed as functions 

of x and u, instead of x and u. Using (3.5) 
l i 
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x(t) 


= <Kt)x. + r (t)u. 

ill 


iT j£ T < (i + 1)T , 


therefore, 


(i+l)T 

J ( x T Ax + 
IT 


u Bu)dt 



[xTuT] 

i i 


L A 21 



where 


11 



12 


o 

T 


0 T (t )A$ (t )dx 


[r ooat^Ct) + B]dT 

3> T ('r)Ar i <-r)dT 


(3.14) 


(3.15) 


A 


12 



d_, Numerical algorithms for calculating the matrices 


A ll* A 12’ A 22 * 


There are various methods for calculating the transition matrix 
$(t) and its first integral (t ) • Most of them are included in 
computer libraries. Efficient and simple algorithms for calculating 
the matrices Q^, A 1X > A 21* A 22’ however > are not found * A highly 

complicated method for evaluating the expressions for A^, A \2 9 anc * ^22 
is described in Reference AS-1. This method requires about 250 FORTRAN 
statements. We will describe a relatively simple method for solving 
(3,15) without a numerical integration, which involves about 80 FORTRAN 
statements. This solution is accomplished in two steps: first, we 

transform all the A f s to a simpler form; second, we subdivide the 
sampling interval T into 2^ subintervals Z£T[£T = (T)/(2 k )] # Then 
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we show that if the various A's are known for [A = A(zXT)], then 

A(T) is obtained in k recursive computations. We will explain 
exactly how it is done and how k is chosen. 

The A's are transformed first because the integrands 

in the integral expressions of these matrices are highly complicated; 

e.g, the complete expression for A is: 

^ 2 l 


21 


T 

= ^ |[j •(T')G 1 d*] [A { (t)GcItJ 


+ B>di , 


(3.16) 


In order to simplify the calculation of A ^ and A 2l' We reformulate 

the cost function of (3.11) by augmenting the state vector x with in- 
put u. The cost function J will be 


J 



uj] * tT (t) 


— 

A 

0 

( t) 

X . 
1 

0 

B 



L J 




(3.17) 


where (x) is defined as 



®(t) 

r x (t) 


'f q x 



- 

= exp 



0 

I 


0 0 


(3.18) 


The expressions for A , A and A 

1 1 1 z ^ z 


are given by evaluating 


i 



A 0 

( x)dT = 

A 11 A 12 


0 B 


A A 




21 22 

_ — 


(3.19) 


By using this formulation we have transformed the problem to one of 
evaluating an integral of the type 
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(3.20) 


Q d = j $(t)Q 0 (T) T dx. 
o 

It is not necessary to calculate integral (3.20) by a numerical integra- 
tion. As in Johnson [JOH— 1], we will use the property of the trans- 
ition matrix 0>(T). 


<I>(T + t) = ®(T)«(x) = 0 >(t)$(T) . (3.21) 


T will be subdivided into 2 parts AT. If Q d (AT) and $(AT) are 
known, then Q(2AT) and 5>(2AT) are given by 


<t( 2 At) 

Q (2AT) 
d 


= 0(At) 1 I i (At) 


2AT 


= ^ <&(x)Q$ T (?)dT 


AT at 

= J $(t)Q0(t) T +J 0(AT + t)q$ T (At + x)dx 


(3.22) 


(3.23) 


Q d (2AT) = Q d (AT) + 4(AT)Q a (AT)$ T (AT) . ( 3 . 24 ) 

k 

Recall that = T/2 , $(T) and Q^OT) are obtained by k recursive 

computations of (3.22) and (3.24), which are relatively simple expressions. 

In order to calculate initial values of O(ZST) and Q^Z^T), 
a constant is selected (say, k = 3) and the following approximations 
are evaluated: 


Q d< AT >i 

<J d (AT > 2 


= qat 
= qat 
®(At) 1 

<&(at) 2 


at 2 t at 2 _t at 1 

+ FQ — + QF — + FQF — 


= I + FAT 


= I + FAT + F 


at 


(3.25) 


(3.26) 
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is smaller than 


If every term of |Q di - Q d2 | and - $ 2 | 

some predetermined number e, then Q^(AT) 2 and ^(Ar)^ , are the initial 
values. If not, k is increased by 1. Throughout our computation 
(short period example), k was never larger than 4. Therefore the 
sampling interval T was divided into 16 parts AT. The predeter- 
mined number e depends on the engineering judgment of the user who 
will have the time constants of the system which interests him. 


B. OPTIMAL DISCRETE REGULATOR 

In the previous section we described how a continuous system, 
control via a ZOH, can be formulated as a discrete system on the 
sampling points. In this section we will develop the theory of an 
optimal regulator of a discrete system which minimizes a quadratic cost 
function. Using the discretization procedure, the theory and results 
of the discrete regulator can be applied directly to a continuous system 
controlled by a digital computer. 

B-l The General Formulation of the Discrete Regulator Problem 

In the last section, it was shown that a continuous system controlled 
by a ZOH with an associated continuous cost function J can be reformu- 
lated into a discrete form: 


X . „ 

1+1 

= Ox. + 
1 

N-l 

r i u i 

fA, „ 

A, _ 

J 

= T 
2 ^ 

r T T i 
[x i u i ] 

11 

12 


1=0 


_ A 21 

A 22 


x . 

l 


(3.29) 


System (3.29) can be transformed into a simpler form. The cost function 
J will be rewritten as 
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= 3 £ [ (u I + x I A 12 A 22 ) A 22 (u i + A 22 A 21 X i ) 

1=0 L 

i i^\l “ A 12 A 22 A 21^ X i] 


+ x. 


(3.30) 


or 


where 


V 


N-l 


j = I y; x^ax. +u t <b < u. 

2 i l 11 


a 

SB 


1=0 


A 11 ' A 12 A 22 A 21 


= A 


22 


U i + A 22 A 21 X i • 


(3.31) 


(3.32) 


Using the last expression of (3.32), (3.29) can be reformulated as 


x. , = (a> - r. a 1 a ) . + r. ii. 

1+1 1 22 21 l 1 l 


1 N-l 

J = - y x ax. + li.BU. 

2 “ i i i i 

i =o 


(3.33) 


If a linear, full state feedback is determined for the system (3.33), 
i.e. , 

II. = Cx. (3.34) 

i i 

then the control for system (3.29) will be 

U i - < C - 4 22 A 21 )x i * < 3 - 35 > 

This transformation is necessary because the theory of the optimal dis- 
crete regulator is solved for the system (3.33), while a discretization 
of the continuous system yields cost function (3.29). 
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For the same reason we can transform a continuous system with an 
FOH to an equivalent form of (3.33). Using our new formulation of FOH, 
(3.2) yields 


x 


i+1 


$x. + [T, 

1 l 




(3.36) 


where I\ is 
k 


B-2 The Solution of the Optimal Regulator 

It was shown that most of the control configurations which interest 
us can be rewritten in the simple form of (3.33) repeated here: 


c 

r = \ $(T)GTdT . 


(3.37) 


*i + i = to i + 


N 


= — Y' x T Ax + u T Bu. . 
2 i i i i 


1-0 


The optimal linear controller is a control law 

u. = C (i )x. (3*38) 

l l 

that minimizes the cost function J, for any initial conditions* If 
N increases to infinity and a steady state is reached, then C(i) = 

C = constant, and the controller is called a regulator. 

The solution of the optimal linear controller was given by Kalman 
CKA-1] who used the dynamic programming approach* We will use Bryson’s 
approach [BR-1] which solves the system (3,33) via the calculus of 
variation* 

For a finite N, the last control is meaningless. It will 

influence only the state x which does not interest us* Therefore 
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the cost function J obtains the form 


* - + 1 1 - “i B “i) • 


(3.39) 


In the minimization procedure used in the calculus of variation, we will 

augment the cost function J by the constraints multiplied by a Lagrange 

T 

undetermined multiplier \ (vector). The constraints are the equations 
of motion for (3.33). The augmented cost function J is 


J = 3 4% ’ a n x n + E (}1 1- Vi 1 + 


(3.40) 


where is defined as the Hamiltonian sequence: 


- \ x I Ax i + i + ^+i (tx i + ru i> 


(3.41) 


Using the methods of the calculus of variation, the condition for 
a stationary value of J is that dj is zero for arbitrary du^: 




+ *r u dx + 1 u du • 

ox o du o 
o o 


We choose such that 


i = 0, . . N-l 


(3.42) 


(3.43) 


T T 

Vn - ■ ° 


(3.44) 


For an extremum: 
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0 


(3.45) 



yields 



T T 

u!b + r = 

i i+i 


(3.46) 


-1 T 

U. = -B r X. . . 
1 1+1 


Combining (3,33)' and (3.47), we obtain 


-1 T 

x. , = $x. — TB r X. , 

l+l l l+l 


and from (3,42) and (2.40), we obtain 


X ± = $ \ +1 + Ax i 


(3.47) 


(3.48) 


(3.49) 


Equations (3.48) and (3.49), called the "Euler- Lagrange difference 
equations" formulated in state space notation are: 


- - 

X 


o + rB -1 r T <&” T A 

-1 T -T 

-rB r ® 


X 

\ 

i+1 

-g>“ t a 

.-t 

$ 


\ 

_ — 


i 


(3.50) 


This is a two— point boundary value problem: is given at i = 0. 

From (3.43) we get the boundary condition for i = N: 


X N = Ax n . (3.51) 

The solution to this problem was accomplished by Bryson using the 
"sweep method" [BR-1] . The sweep method assumes a solution for \ 
of the form: 

= S x . (3.52) 
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This solution leads to a matrix Riccati difference equation in S : 


S . 
J 


T -1 Tv 

<& (S + TB T ) 0 + A 


j = N-l, o 


S N = A ' 


(3.53) 


Determining S. from the backward recursive relations (3.53), and using 

(3.57) and (3.42), the optimal control is expressed as a linear 

combination of the state x ± . If certain conditions are satisfied as 

N increases, S. reaches a steady state S , and the controller is 
J ss 

induced to a regulator: 


u = -b” 1 T$" T (S - A)x . (3.54) 

i ' ss i 

which is obtained by combining; (3*47), (3*49), and (3.52). In the steady 
state, the matrix Riccati difference equation is reduced to a second- 
order matrix algebraic equation 

S = $ T (S -1 + rB -1 r T ) $ + A. (3.55) 

ss ss 

During the last two decades, a considerable effort has been made to find 

an efficient solution of the Riccati equation and the steady state 

matrix equation. The usual method of solution for (3.55) is a recursive 

computation of (3.53) until S reaches a steady state S . A com- 

s s 

pletely different approach to solving for S is to use the eigenvector 

s s 

decomposition of the transition matrix (3.50). 

3. Solution of S by Eigenvector Decomposition 

' u "™ ' L "“ 1 “ " " 3 S — - - 

In 1966, Potter [PO-1 ] described a method for the steady state solu- 
tion of the matrix Riccati differential equation by eigenvector decomposi- 
tion. Bryson and Hall [BR-2], using efficient QR algorithm for eigen- 
vector calculation, constructed a computer program for linear regulators 
and Kalman filter syntheses. Vaughan [VA— 1 ] extended the Potter method 
tor discrete system control synthesis. The authors solved the eigenvector 
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decomposition problem independently of Vaughan, and also applied it to 
the discrete filter synthesis problem. Our interpretation and our proof 
will be given here. 

The Euler— Lagrange equation (3,50) will be repeated here, but on 
the z-plane: 


zx 


z\ 


-1 T -T - 1 T -T 

o+pb r <t> a -ra r <& 


-T 

-0 A 


-T 


The following theorems will be proved. 


Theorem 1 : If z is an eigenvalue of the system (3.51), the l/z is 

also an eigenvalue. 


Proof : (a) defining E: E = z ; (b) defining a new variable: y = z\ 

and directly using (3.48) and (3.49), the system (3.56) can be 
transformed to an equivalent form: 


-1 T 



$ - Iz -TB r 


X 

T -1 



A $ — Iz 


r 



- - 


0 


(3.57) 


-1 x 

TB T and A are symmetric; therefore, (3.57) may be rewritten 
as 


T T 

Cr x ] 


T 

<t>- ie -pb r 


T -1 
A 4> -IE 


= 0 , 


(3.58) 


Systems (3.57) and (3.58) have the same transition matrix. There- 
fore, if z is the solution of the characteristic equation of 
(3.57), then E is a solution also. 


-2 7- 



Conclusion : The eigenvalues of the Euler- Lagrange equations (3.56) are 

reflected symmetrically across the unit circle on the z-plane. See 
Fig. Ill— 2. 



FIG. 1 1 1-2 ROOTS LOCATION OF EULER- LA GRANGE (E-L) 

EQUATIONS 


Definition: 



= f T z V’ 


(3.59) 


Matrices T and T are the eigenvectors of the E-L equations, 
z t 

(3.56), associated with z and E respectively. 

Before formulating and proving Theorem 2, a well-known result from 
linear systems theory will be presented [e.g., KW-1]. 

A homogeneous, linear time-invariant discrete system 

x. _ = $x. (3.60) 

i+l l 

with the initial conditions x , has the solution 
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x. 

1 



(3.61) 


The solution x. can be expressed in terms of the individual eigen- 
vector modes as follows: the initial condition is resolved along 

the modes of the eigenvectors by the transformation : 

g = T -1 x (3.62) 

• o r O 


where T is the matrix of the eigenvectors of 0 
r 

T r “ • < * > j t * * • t ] * 


(3.63) 


The solution x. is a linear combination of the particular excited 
modes , i.e. , 


x . 

x 



Using this result, we can formulate Theorem 2. 


(3.64) 


Theorem 2 : The steady state solution S gs at the matrix Riccati differ- 

ence equation (repeated from Eq. 3.55) 

S = $ T (zS -1 + rB -1 r T ) -1 3> + A (3.65) 


is 



(3.66) 


Proof: The homogeneous solution of the E— L equations (3.56) is 


x 


X 


i 




(3.67) 
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where n = the order of the system, 

( ^ °k 

< — constants expressing boundary conditions 


= eigenvalues of (3,51) 


x (k) 
z 


\ (k) 

z 


*E (k) 


V k> 


are the eigenvectors of (3.56) corresponding to z^ and E fc . 
tion (3.67) may be formulated in a matrix notation 

r i r 1 r 1 


(1) (n) 

X X 

z z 

(4) (n) 

z z 


(1) (n) 

E E 


^i z i 


l 

s on n 


defining Z as 


z i 0 


0 z n 


Equation (3.68) can be further reduced to 


x i = V So + V \ 

\ = A z Z l + V N_i \ • 

As i increases, the stable modes multiplied by the vector | 
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attenuate so that equations (3,70) become 


„ „N-i - 

X E Z \ 


vr \ 


Solving (3.71), we get 


X. = AjCV . 
1 E E i 


(3.71) 


(3.72) 


But \ = S x was the assumed solution of the matrix Riccati 

i i i 

difference equation (3.53). Therefore (repeat of Eq. 3.66) 

S s5 = Ve 1 - 

This concludes the proof of (3.66). 


Having these results, the optimal feedback control of the linear 
discrete regulator is 

u = -B _ 1 r<D~ T (A x ” 1 - A)x. = Cx. . (3.73) 

i EE i l 

This solution requires nonsingularity of the transition matrix 
Made discrete, the linear continuous system always has the property 
that |o| ^0. This stems from the fact that a continuous linear 
system has a unique solution for a given initial conditions [e.g., 

KW-lJ; however, this ‘ property ( |o | ^ 0) is not obvious for a pure 
discrete system. If, in a pure discrete system $ is singular, it 
indicates that some of the states (or the modes) can be expressed as 
a linear combination of the remaining states (or modes). Hence, if a 
state variable feedback is required, the singular system representation 
can be reduced in dimension until regularity is achieved and the feed— 
back matrix C can be calculated. 
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Two more results will be given: 


a) The stable eigenvalues of the E-l equations (3.56) are identical 
with the eigenvalues of the closed loop optimal system. This useful 
property can be proven by analogy to continuous systems [e.g., BR-2]. 

b) The expression for S , Eq. (3.66), repeated 

ss 

m 

S s s = Ve 1 

is independent of any rearrangement of the individual eigenvectors in 

- - 

-V 


The proof is obvious if (3.66) is rewritten as 


S 

ss 


A 


E 



(3.74) 


Assuming 
tion of 
relative 


S is fixed f then any column j of 
s s 

S and column j of A_. Its values 
ss E 

position with respect to other columns. 


are 


is a linear combina- 
independent of its 


The computer program, DISC, listed and explained in Appendix A, 
is based on these results. 


£. THE OPTIMAL LINEAR DISCRETE FILTER 
C-l The Measurement Timing 

The purpose of the filter is to reconstruct the states which are 
not measured, and to minimize the process and measurement noise influ- 
ence. For the physical system, 


1+1 

- $x. + r. u, + r w. 
1 11 2 1 

w ± ■+ N(0, Q d ) 

y i 

= Hx. + v. 
l i 

Vj . -♦ N(0,R) . 
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The optimal steady state filter is given by 


i+l 

r= Ox. + T.U. 

(3.76) 

1 11 



II 

X 1 

+ 

* 

s. 

< 

1 

X 1 

w 

(3.77) 

i+l 

i+l i+l i+l 

T -1 


K 

= PH R 

-1 T v -T 

(3.78) 

P 




where M is the error covariance matrix before measurements, and P 
is the error covariance matrix after the measurements y^ + ^ were done. 
How M is computed will be explained later. Relation (3.77) is the 
one which is used by Bryson and Ho [BR— 1 ], and throughout this work. 

Here, we assume a zero computation time between the measurements y^ ^ 
and the output of the filter x . Borow [B -1] shows that the delay 
for the short period mode calculation is of the order of one millisecond. 
The other possibility is to use the approach of Kwakernaak [KW-1]: 


i+1 


x. + 
i+l 


K(y. - H$.) 
'i l 


(3.79) 


These two approaches could be summarized on the time axis in Fig. 
Ill— 3. 


MEASUREMENT OUTPUT MEASUREMENT OUTPUT 

\ \ f 



iT (i+l)T 


+- 

t 


FIG. II I— 3 


MEASUREMENT AND SAMPLING INTERVAL 



In (3.77) we assume e *+ 0. Both of the formulations have 
advantages and disadvantages. The main feature is that (3.77) is an 
approximation only as e is finite. On the other hand, (3.69) uses 
’obsolete* information. The compromise will be to use a third method; 
for example, the one shown in Fig. I I 1-4. 


M EASUREMENT OUTPUT 

K v 



j 

t 


FIG. 1 1 1-4 INTERSAMPLE MEASUREMENTS 


6 q is the fixed time interval, longer than the maximum computation time 

needed to generate the states _ from y. „ measurements. As is shown 

l+l 'i+l 

by Kwakernaak et al [KW-1], the calculation of the optimal filter is 
highly complicated. We will use the first approach which assumes that 
the delay time between the measurements and the output of the filter £ 
is much smaller than the sampling. From an inspection of (3.76) and 
(3.77), we can see that most of the updating of the filter can be done 
before the measurements y. are received* Actually, the only calcula- 

X + X 

tion which has to be done during the delay interval is to multiply 
K by and to ad d this quantity to the rest of (3.77). 

C-2 Calculation of the Steady State Optimal Filter by Eigenvector 
Decomposition 

The steady state optimal (Kalman) filter is an observer which mini- 
mizes the steady state error covariance matrix P. The recursive equa- 
tion for M obtained from the minimization process has the same struc- 
ture as those for S for the optimal control problem. The equations 
are identical if we consider the following equivalence first recognized 
by Kalman [kA-1]: 
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Control 



Filter 

<3> 

— 


_t 

$ 

H 

S' 


i 

V" 

x _ 


r ! 

h t 

V." 

s 


A 

Q d 

A 

11 


B 



R 


Thus we will replace the matrices of the optimal control problem by the 
matrices of the filter in E-L equation (3.56) which yields the matrix 


K: 


H = 


T T —1 —1 T 

0 + H R m r 2 Q d r 2 


-1 T 
r 2 Q d r 2 


T -1 
-H B 


-1 

4 * 10 



(3, 80) 


Making direct use of the results from the optimal regulator calculation, 
the error covariance matrix M is given by 


M 


x e a e 


where X_ and A are defined from 
E E 



(3.81) 


(3.82) 


T is a matrix of the eigenvectors of the 2n X 2n matrix X , 

X 

z 

A 

z 
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are the eigenvectors of H corresponding to the stable eigenvalues of 
K, and 



are the eigenvectors corresponding to the eigenvalues of X located 
outside the unit circle. 

Similar to the optimal control problem, the stable eigenvalues of 
K (inside the unit circle) are identical with the eigenvalues of the 
observer error system defined from (3,76) and (3,77) as 

x £ S - x (3.83) 


X 


i+1 


(0 - KH#)x. + (K Hr - r )w. + Kv. • 
i 2 2 i x+1 


(3.84) 


D, ALGORITHM FOR AN EVALUATION OF THE STEADY STATE 
RESPONSE TO AN EXTERNAL NOISE 

A stable discrete system disturbed by an external noise with 

a covariance matrix Q t reaches a steady state. The average behavior 

d 

of its states is characterized by a covariance matrix X which is the 
solution of Bryson and Ho [BR-1] : 

X = $X4> T + Q d . (3.85) 

Following Bryson, the average behavior of an optimally controlled dis- 
crete system, with an external noise disturbance and measurement noise, 
is characterized by the state covariance matrix X. This is the solu- 
tion of 

x - m = ($ + rc)(x - p)(r + rc) T (3. 86) 

where C is the optimal gain (3.73), M and P are the observer error 
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and covariance matrices defined in (3,78)* Equation (3,87) will be re 
written as 


X - P 


(o + rc)(x - p)M> + rc) 


T 


4 M - P . 


(3,87) 


Equation (3,87) is now in the form of (3,85), Equation (3,85) is essen- 
tially a linear equation of X, New algorithms for solving (3,85) 
appear frequently in the numerical method literature [e,g, , BER— 1], V^e 
will present a new numerical solution of (3,85) which utilizes the al- 
gorithm for eigenvector decomposition. 

Claim 1 : X - (3,88) 


where 



is a matrix of eigenvectors of the 2n X 2n system of H: 


/ 

H 




(3,89) 
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is the eigenvectors submatrix corresponding to the eigenvalues of H 
outside the unit circle. 

' -1 

Claim 2: If z, is an eigenvalue of H, then z,_ = E, is also an 

■ It K rt 

/ 

eigenvalue of H. 


Proofs : Claim 2 is obvious from inspection. To prove Claim 1, we 

let 




where 


D 


"V 

D 




> 1 . 


(3.90) 


(3.91) 


Using (3.90) and (3.91) 


= Ve 

V" Tx e + m e ‘ Ve 

yields 

Q d + *Ve V = V * 1 ■ 


(3.92) 


(3.93) 


Define 
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X 


(3.94) 


A E X E ' 

Equation (3.93) is identical to (3,85), thus proving Claim 1. 

Using this result, the average behavior X of the optimally con- 
trolled discrete system, including process noise and measurement noise, 
is given as the solution of (3,81), which is 

X = AgX” 1 + P , (3.95) 


where 




a e 


is a matrix of eigenvectors corresponding to 2n X 2n matrix H; 


H = 


(5> + Tc) 


-T 


(m - p)(o + rc) $ + rc j 


(3,96) 
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E. SUMMARY 


1* A linear continuous system, controlled by a discrete controller 
via a zero order hold (ZOH) or a first order hold (FOH) can be formulated 
as simple linear discrete systems (3.97) # We have shown that if a digital 
computer generates the control sequence, the information needed for the 
FOH can be generated simultaneously, and not by extrapolating the last 
two control inputs* 

If a continuous cost function is assumed and a continuous white 
noise acts on the system, the whole controlled system can be represented 
by equations of the form 


X i+1 


i + IJ_ U ± + r 2 W i w i -* N (°> Q d ) 

N-l _ 


J = 


£ x i ta i + u i Bu i 

i=o 


y. = Hx. + v. 

1 IX 


v ± -» N (°^ r) 


(3.97) 


Converting to the formulation (3*97) enables us to calculate the optimal 
^regulator and the optimal steady state observer using discrete algorithms. 

2. We have developed an eigenvector decomposition computer program 
for solving the discrete steady state matrix Riccati differential equa- 
tion. It is highly efficient compared to recursive methods, due to the 
use of the QR algorithm [FR-l], which finds eigenvalues and eigenvectors 
of a matrix with widely dispersed eigenvalues very rapidly and accurately. 
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Our eigenvector decomposition program calculates the optimal 
regulator gains and the optimal steady state filter gains, and as a by- 
product, it provides the closed loop eigenvalues and the observer error 
poles [ App # A]. 

3. Calculation of the root mean square response of the closed loop 
system involves the solution of a linear matrix equation of the form 

X = 0X0 T 4 Q d . (3.98) 

A new algorithm for solving (3.98) was found. This algorithm is based 
on eigenvector decomposition of a system associated with (3.98). We 
do not claim that this algorithm is more efficient than algorithms 
already in the literature* Its main advantage is that it uses the 
existing program for eigenvector decomposition of a matrix. 


- 41 - 



IV. SAMPLING TIME AND TIME RESPONSE 


The time response of a controlled system is determined by its 
closed loop dynamics and by the input signals, A system has a proper 
time response if for a specific input , some combination of the states 
follows approximately a predetermined pattern. In the case of a pilots 
commands, the input is a continuous signal, e.g., an output of a poten- 
tiometer. This continuous signal has to be sampled and fed to the 
digital computer which further processes it. The questions now are: 
how well does the digital processor interpret the sampled signal, 
and how well do the states follow the desired pattern? 

The theoretical basis for the reconstruction of a continuous signal 
from sampled data is given in the sampling theorem developed by Shannon [e.g., 
RA-1], , The sampling theorem states that in order to reconstruct an unknown 
continuous signal from samples of that signal, one must use a sample 
rate which is twice as high as the highest frequency contained in the 
unknown signal. This theorem is not directly applicable to a reconstruc- 
tion of an input signal in a real time digital system. This is true for 
the following reasons: (a) for causal systems, the reconstructed signal 

has a phase shift with respect to the input signal. The theorem states 
that it is possible to reconstruct the signal but it doesn’t say that 
it will be done in the same time. Actually, any reconstruction scheme 
has to accumulate a minimal amount of data points in order to start the 
reconstruction of the sampled signal. (b) For aircraft applications, the 
proper time response is* formulated as a response to the pilot's step in- 
put. However, from a theoretical point of view, a step input has an infi- 
nitely large frequency spectrum. (c) In the feedback path, the sampling 
theorem doesn't hold since we have knowledge of the plant. Therefore, in- 
steady of applying the sampling theorem, we will investigate the time re- 
sponse of a realistic example using a simulation scheme. 

The key question of how to get a proper time response can be reform- 
ulated as a problem of what the input signal should be in order to 

*** m £Imm 
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obtain some predetermined pattern of the states. Extensive research 
dealing with this complicated problem is currently being done by Holley 
and Bryson [HOL-1]. Their approach, called ,f a nonzero set point regulator 
design", consists of an analytical basis for the time response synthesis, 
and can be extended to discrete control systems. Meanwhile, we will use 
the time response specification given by Borow et al [BO-1] and Sutton 
et al [SU-1], which is based on the accumulated experience of the U.S. 
Navy. A complete description of this specification will be given in 
Section IV-B. 

The outline of this chapter is as follows. In Section IV-A we 
will define the example, "F-H", which will be used in this work. It 
includes the short period mode, a bending mode, and wind gusts. 

In Section IV-B we will describe the various design objectives, 
including the time response. 

In Section IV-C the inner control loop design will be described. 

The classical approach will only be outlined, but the optimal discrete 
design will be explained in detail. 

In Section IV-D the main objectives of this chapter will be inves- 
tigated. It will be shown that, for the F-H short period example, a 
proper time response imposes a definite limit on the sampling interval. 


A. EXAMPLE DEFINITION 


In order to be able to illustrate the various factors which influ- 
ence the sampling rate selection, our investigation into aircraft discrete 
controls should be related to a definite and relevant technical system. 

We chose a controlled short period mode configuration of a hypothetical 
type aircraft, which we term F-Hi, This type of aircraft is described 
by Borow [BO-l] and Sutton [SU-l] as the future aircraft of the U.S. Navy. 
Our choice was made for the following reasons: 

1. The time constants and the time responses in the longitud- 
inal mode are short and they constrain the sampling rate. 

2. The influence of bending modes on the stability and sensitivity 
of the control loop must be taken into account, A considerable 
effort must be made to include these effects in the control de- 
sign. 
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3. Only a military aircraft flies in such difficult flight condi- 
tions as Mach 1.2 at ground level. The strongest wind gust 
amplitudes are found at zero altitude. A high velocity flight 
through such gusts generates a short correlation time of the 
external disturbance, which makes the gust alleviation more 
difficult. During landing conditions (Mach 0.19), this type 
of aircraft has an extremely slow time response (oj = 0.5 rad/sec). 
Therefore a strong bandwidth augmentation is necessary. It will 
be shown that this behavior imposes serious limitations on the 
sampling rate. 


The geometry of this aircraft as a flexible body is described in Fig. 
IV-la and Fig. IV-lb . References BO-1, ED-1, and SU-1 use this model. 

The equation of motion, based on body axes and dimensional stability 
derivatives are 

q = (M )q + (M )a„ + (M.)a + (M. )6 
q a T a T 6 e 


= Q + ~ a T + 


w 


fii 0 

y W W 


U T u 
o w o 


6 - * " — - 71 

e x U 8 
w o 


w 


w /§r a 

g , V W w 
-r 

T T 

W W 


Tj 

'n 


(4.1) 


These are the short period mode equations; they include a vertical gust, 
modeled as a first order Gauss-Markov process. 



(4.2) 


is defined as the total angle of attack. The reason for preferring the 
state variable instead of a (as used by Borow [BO-1] and Sutton 

[SU-1]) is that 0^, is an important variable in the gust alleviation 
study. The vertical acceleration is directly related to 

The bending mode is modeled as a second order system driven by the 
elevator input and the rate of the angle of attack — a. The bending mode 
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is slightly damped (£ = 0.01). The bending equations are 

b 


X 3 = W b X 4 


x . = -CO, X - 2£, to X + co, k 2 5 + k o co Z a • 

4 b 3 b b 4 bl5e 2 b a 


(4.3) 


The actual acceleration and rotation of any point on the body axis 
to bending, is 



X 1 dUe 


(4.4) 


where k^, k^, are quantities depending on the airplane shape and mass 
distribut ion, and z^(x) is the first mode shape, .The wind gust influ- 
ences the bending mode primarily through 6^ (see Eqs. 4.1 and 4.2). 
Sutton [SIT-1 ] and Blakelock [BL-l] disregard this influence* 


An accelerometer and a rate gyro are the measuring instruments. 
The measurement signals are combined from the rigid body motion, the 
bending mode motion, and additive white noises: 


Q 


T 


n 


T 



v 

q 


U ( q - (5) + X q -h tt— x A + v < 
o ^ oj 4 n 

b 


(3.5) 


where v = N(0.r ) . v = N(o,r ). For a numerical example, we chose two 
q 7 q 7 n 7 n 7 

difficult flight conditions of the F-H aircraft [B0-1 ; SU-l]; 

a) Flight condition No. 1 is a flight at zero altitude and Mach 
0.19. As mentioned in the beginning of this section, a con- 
siderable increase in the bandwidth (acceleration feedback) is 
necessary in order to improve the time response. It will be 
shown later that this flight condition limits the selection of 
the sampling rate. 
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b) Flight condition No. 8 is a flight at zero altitude and Mach 
1.2. High velocity flight at zero altitude generates strong 
wind gusts. The maximum available gust alleviation is neces- 
sary for several reasons, including improved aiming and mini- 
mizing the fatigue of the pilot. 

The corresponding numerical data are summarized in Table IV-1. 


Table IV-1 

NUMERICAL DATA OF THE F-H SHORT PERIOD ( SP) MODE 



Condition No. 1 

Condition No. 8 

U ( ft/ sec) 
o ' 

212.0 

1340 

M ( l/sec) 

-0.44 

-1.91 

<1 



M (l/ft sec) 
w 

-0.0017 

-0.13 

M &e ( l/sec 2 ) 

-1.23 

-69.1 

M-U/ft) 

-0.62X10* 3 

0,5 2x1 3*” 2 

Z w ( l/sec) 

-0.57 

-4.03 

Z &e (f l/sec 2 ) 

-13.2 

-399.0 

t (sec) 

Ml 

3.0 

0.5 

CT (ft/sec) 
w' ' 

12.0 

12.0 


o 

• 

o 

0.01 

ua^( rad/sec) 

25.0 

25.0 

6 ( rad) 

0.4 

0.4 

z^ of accelerometer 

0.07 

0.07 

of rate gyro 

-0.005 

-0.005 

k i 

4.0 

4.0 

k 2 

0.06 

0.06 




l (1) 

14.0 

14.0 

a 



SP freq (rad/sec) 

0.60 

13.2 

period of SP (sec) 

10.5 

0.476 
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A- 1 Open Loop Coupling Between the Bending Mode and the Rigid Body 

Motion 

Borow et al, Blakelock, and Sutton et al [ BO-1, BL-1, SU-l] dis- 
regard any direct influence of the bending mode on the rigid body motion. 
However, in flight condition No. 8, the supersonic velocity in a high 
density atmosphere generates a considerable restoring moment (M^) . Conse- 
quently, the short period mode oscillation is high (~ 13.6 rad/sec) and 
has the same order of magnitude as the bending mode frequency (25 rad/sec). 

To estimate the influence of the bending mode on the rigid body 
motion, we will assume a simplified configuration which generates a moment 
around the center of mass, see Fig. IV-2. 



FIG. IV-2 MOMENT AROUND THE CENTER OF MASS AS GENERATED 
BY THE BENDING MODE, 


The force F^ and the corresponding moment around the center of mass are 
proportional to z'(x). , We may consider the rotation of the tail as an 
additional angle of attack of the elevator. In supersonic aircrafts, the 
whole elevator surface is moving and thus the moment caused by F is 

approximately equal to -Mg (z'x )/(u> ) where <z ! x )/(</) is the 

6 G u D G 3 u 

additional angle of the elevator due to the bending rotation. 

We will use another assumption: the sum of the moments generated 

by the deflection of other parts of the body ( including F n ) is equal or 
lower than 50$ of the moment generated by the bending rotation of the 
elevator. 
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By using these approximations, the additional term of the first 
equation of (4.l) will be 


“ k b M 5 


Z e X 3 


where 0 ^ k £ 1.5 . 
b 

The poles of the coupled system, combined from the rigid body and 
the bending mode are summarized in Table IV-2. 


Table IV- 2 

THE OPEN LOOP POLES OF THE COUPLED SYSTEM 



Rigid Body 

Bending Mode 

k b = ° 

s = -2.5 ± j 13. 6 

s = -0.5 ± j25.0 

k = 1 

s = -0.9 ± J13.5 

s = -1.9 ± j24.8 

b 


\ = 1 - 5 

s = -0.3 ± j 13.4 

s = -2.4 ± J24.7 


From Table IV-2 we see that for a detailed design, the coupling be 
tween the short period mode and the bending mode cannot be neglected. 
Although the stability augmentation system stabilizes the system, the 
open loop configuration is only marginally stable. 


B. DESIGN OBJECTIVES 

There are several different objectives required of an inner loop 
longitudinal autopilot. The most important may be summarized as follows 
( l) A proper time response to various inputs; (2) A proper dynamic be- 
havior based on the pilots experience; (3) Wind gust alleviation. 
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The crucial constraint on the actual design is the condition that 
for one definite flight condition, there should be only one control con- 
figuration which will meet all the objectives. 

For the conventional aircraft where the elevator is behind the 
center of the mass, the objectives 1, 2, and 3 are conflicting; thus, 
a suitable compromise will be suggested in Section IV-3, 

The various design objectives will now be explained in more detail. 
Objective 1: The proper time response requirement varies for different 

aircrafts. Essentially, the requirement is for a fast response, but one 
which is not too sensitive and has a sufficient stability margin. The 
exact formulation will be given with Objective 2, 

Objective 2: A desirable dynamic behavior, expressed as a location 

of the augmented short period poles, is described in Kolk [KO-1], The 
shaded area in Fig, IV- 3 is the location of poles which is preferred by 
pilots , 



FIG. IV- 3 SHORT PERIOD POLES LOCATION; PREFERRED 
BY MOST PILOTS. 
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Combining Objectives 1 and 2, we may say that a high performance 
military aircraft is expected to have a fast response in order to be 
able to execute properly different maneuvers. It must also dampen ade- 
quately for precise weapons release. An attempt to give a quantitative 
formulation to these requirements was made by the Navy, as described by 
Borow [BO-1] and Sutton [SU-1] under the name TRP (time response parameter). 

The TRP is defined by the quantities in Fig, IV-4 and the relation, 


TRP 



+ 0.08 (A -1) + 0.05 (t H 

q a nz 


0.7) 


+ 0. 3 (A - 0.3) 
nz 


(4.6) 


where A is the pitch rate overshoot, and A ' is the vertical accel- 
q 7 nz 

eration overshoot. The objective is to keep TRP < 0.25. All parenthe- 
sized terms, if negative, are assumed to be zero for the TRP calculation. 



FIG. IV- 4 THE TIME RESPONSE PARAMETER (TRP) 

DESCRIPTION 
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We will interpret the meaning of the different terms in the TRP 
specification: TRP is zero (the 'best 1 response) if: (a) the overshoot 

of q is not larger than 1; (b) the time to reach half final value of 

the vertical acceleration is less than 0.7s; (c ) the maximum overshoot 
of the vertical acceleration is less than 0.3; (d) the undershoot (the 

non-minimum phase property) is shorter than 0.2 sec; (e) the system is 
well damped. Therefore, (t c )^ shouid be as large as possible, but not 
lower than (t_) = 4(td). 

c q 

Generally speaking, the TRP specification indicates that the vertical 
acceleration time constant should be less than ~ 1 sec* 

Objective 3. Gust alleviation is achieved by a proper control* 

The optimal controller, which minimizes a quadratic cost function 
for an impulse response, minimizes also the mean response to white noise 
disturbance (Parseval Theorem [KW-l])* Further alleviation can be 
achieved if the external disturbance has a finite correlation time and 
if the average wind gust behavior can be predicted* In Chapter V, a 
detailed description of relations between gust alleviation and the samp- 
ling time will be given* 


C. INNER CONTROL LOOP DESIGN 

The main purpose of this section is to point out the assumptions 
and the basic differences between a classical design and the modern 
control approach. Only a brief description will be given. 

C-l Summary of Classical Control Design 

The classical control design of the inner loop can be subdivided 
into two basic methods: (a) A continuous design on the s-plane and 

discretization; ( b) A design in the z or w-plane. Both of these 
methods use the transfer function approach- This classical design is 
summarized in Fig. IV-5. 
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DIGITAL COMPUTER 

Fig. IV- 5 CLASSICAL DESIGN OF THE CONTROLLED SHORT 

PERIOD MODE 


The first method consists of three basic steps: 

1. A continuous design of the inner loop based on previous exper- 
ience. 

2 . Discretization of the compensation, using digital filter approx 
imations ; 

3. Modification of the different gains and parameters, in order 
to diminish the unwanted properties of the sampling (e.g. 
time lags). 

Step 3 can be omitted if the sampling rate is fast enough for the digital 
filter approximations to cause negligible error. It often results in 
5 to 10 times faster than the bending modes. 

Different elements, mechanized on the digital computer, are 
designed on the s-plane by classical techniques. Discretization is 
achieved by the Tustin Transform [BO-1, ED-1, and OS-l] or, in rare 
cases, by z-transform [SU-l]. This design, which works quite well with 
analog elements, has the following disadvantages: (a) Gust alleviation 
can only be achieved by a trial and error computation and simulation. 

(b) The sensitivity behavior on the s-plane doesn’t coincide with the 
sensitivity of the discrete controller. (c) Approximations in digital 
filters lead to undesired properties of the closed loop system (as 
aliasing) . 






The second method consists of two basic steps: 

1. Transformation (discretization) to z- or w-plane. 

2. Classical design of the digital compensator based on previous 
experience. 

The compensation network is usually designed on the w-plane [ ST-1, LE-l] 
because of its similarity to the s-plane (c.f,, Chapter II) . This is 
an exact design and therefore the problem of digital approximations is 
eliminated. 

The main disadvantage of the classical approach is its inability 
to handle multi-input, multi-output design problems. Even Lee [LE-l], 
who is currently developing a w-plane automatic synthesis program, 
says, "The use of a classical approach is limited. Sooner or later, 
we will have to start design in the state space." 

C-2 Optimal Discrete Design 

An optimal discrete synthesis, based on quadratic criteria, includes 
a full state-variable feedback and an optimal steady-state observer. This 
is an exact design; no approximations are made as in the first 
classical design, where the continuous compensation was made discrete. 

The discrete compensator, calculated for a predetermined sampling 
interval, always yields a stable system which minimizes given quadratic 
criteria. However, these properties are only correct if the assumptions 
about the controlled system are perfect; i.e., the system is linear, 
there are no limitations on actuator bandwidths, the designer has a good 
knowledge of the systems parameters, etc. 

Furthermore, for a real system, the designer has to face a whole 
new set of problems: (a) What model to use; (b) How to determine the 
cost function; (c) How to achieve a proper time response; (d) How 
to handle nearly undisturbable states; (e) How to reduce the sensitivity 
to parameter variation; (f) What sampling rates to choose in order to 
obtain a minimum response to external disturbances; (g) How the system 
behaves for longer sampling intervals (roughness). The last three 
problems will be answered in Chapters V, VI, and VII. The first four 
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problems will be partially answered by outlining the optimal discrete 
design of the principal example ♦ The optimal design is based on the 
Separation Theorem [GU-l] which shows that the controller and the ob- 
server can be designed separately. 

a. The controller design . 

The main difficulty in an optimal control design is to find the 
numerical value of the weighting matrix A of the quadratic criteria. 

This problem is further complicated by the nonminimum phase property 
of the F-H short period mode. 

The necessity for a fast response conflicts with the effort 

to alleviate the wind gusts. Therefore, the root square locus approach 

was used [BR-4]. The short period poles are relocated to an acceptable 

region, while in the meantime the gust response is minimized. No weight 

is given to the bending states x^ and x^, since no attempt will be 

made to control the bending mode (during the first design). Similarly, 

no weight will be given to w . 

g 

The root square locus of the short period mode, based on Eq. 
(4.1) is: (a) for a continuous system 

1 + ~b y q (“ s)y q (s) + T y a^“ s ) y a ( ' s ) = 0 ; (4. 9) 


( b) for a discrete system: 

A 


1 + B 2 y q (z ’ 1)y q (z) + T y a (z “ 1)y a (z) 


(4.10) 


where 


y - f 

q 6 e 


a 

y a ” T” * 

e 


(4.11) 


A and A are the weighting terms in the main diagonal of the weighting 


q a 

matrix A. 
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In order to obtain a better engineering insight, the root 
square locus was also traced on the s-plane. Note that the trace on 
the z-plane depends on the sampling interval. 

The root square locus of the F-H short period mode for flight 
condition No. 8 is given in Fig. IV- 7. As can be seen, the minimization 
of a increases the bandwidth, which is undesirable. Therefore, a com- 
bination of A and A was chosen as a nominal design. 

Q Of 

The root square locus for flight condition No. 1 is given in 
Fig. IV- 9. The major problem in this flight condition is the sluggish 
response to the pilot's commands. Therefore, a bandwidth increase was 
necessary. We will analyze the time response in moredetail in Section 
IV-D. 

The nominal pole design for flight condition No. 1 vs the 
sampling interval is traced in Fig. IV-9. The damping is nearly unchanged. 

b. The observer design . 

An observer design, which uses steady state optimal techniques, 
is straightforward if the power spectral density matrices of the noises 
are known. But the nearly undisturbable modes are a problem which needs 
further clarification. A nearly undisturbable mode, in our case, is the 
bending mode, which is primarily excited by the elevator's input [BL-1, 
B0 - 1, ED-l]. In this case, the observer error corresponding to the 
bending mode is virtually undamped. Only Sutton [SU-l] assumes a slight 
excitation of the bending mode by the external disturbance. Consequently, 
this particular mode, the optimal filter yields a very low gain in 
the observer error equations. In order to avoid this situation, an 
artificial noise will be applied to the bending mode state. This is 
done only for computational purposes and it is equivalent to a pole 
placement of the observer error. Another method, which relocates the 
observer's error poles of the undisturbable modes and uses an optimality 
criteria, is currently being developed by Breza and Bryson [BRE-l], This 
method can be extended to a discrete system. 
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FIG. IV-6 POLES LOCATION VS WEIGHTS OF THE F-H SHORT 
PERIOD MODE. Flight Condition No. 8. 



FIG. IV-7 POLES LOCATION VS WEIGHTS OF THE F-H SHORT 
PERIOD MODE. Flight Condition No. 8. 
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FIG. IV-8 POLES LOCATION VS WEIGHTS OF THE F-H SHORT 
PERIOD MODE. Flight Condition No. 1. 



FIG. IV-9 THE POLES OF THE NOMINAL DESIGN OF THE F-H 
SHORT PERIOD MODE vs SAMPLING INTERVAL. 
Flight Condition No. 1. 
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In this chapter we are interested in the time response to the 
pilot 1 s commands. Those commands are fed simultaneously to the plant 
and to the observer. Therefore, during the first seconds, the time re- 
sponse of the plant is unaffected by the observer error poles. In 
Chapter VI, where we deal with the sensitivity problem, the location 
of the observer error poles will be investigated in detail. 


D. THE TIME RESPONSE OF THE F-H SHORT PERIOD MODE EXAMPLE 


One of the objectives of the optimal discrete controller is to mini- 
mize the average pitch rate q and the average total angle of attack 
G^. This is accomplished by determining a proper weight on the states 
in the quadratic cost function. But we are limited in our choice, de- 
pending on: (a) the bandwidth of the actuator ( ~ 45 rad/sec) ; (b) the 

maximum amplitude of the actuator (0.4 rad); ( c) a proper time response. 

The closed loop poles are far below the bandwidth of the actuator 
and the restriction is the time response. We will explain in detail 
how a proper time response is achieved and how the time response is 
related to the closed loop poles and to the sampling interval. 

The link between the pilot and the behavior of the aircraft is 
described in Fig, IV-10. 


PILOT/ 

INPUT 


1 


LOW PASS 
FILTER 


OH 


COMPEN- 

SATOR 


I 


] 1 - 

DIGITA~L COMPUTER 


ZOH 


AIRCRAFT 


FIG. 4.10 THE CONTROL CONFIGURATION OF THE SHORT PERIOD MODE 
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The low pass filter is necessary and serves two purposes: (a) It 

is known from previous experience [BO-1, SU-l] that pilots favor a 
smoothed stick input with a time constant of 0.4 to 0.6 sec. (b) The 
low pass filter helps to reduce the pitch rate overshoot of the controlled 
configuration. 

The only way to visualize the time response to the pilot input is 
to simulate it. Our results are based on a simulation scheme described 
in Appendix B. 

We will now explain in detail how the pilot command is actually 
processed by the system. This detailed analysis is important. 

The time t = 0 will be defined as the instant when the pilot 
executes the stick input (a step function). The information about 

this command will reach the computer within T time (T sampling inter- 
val). In order to be on the safe side, we have to assume a full delay 
T^ = T between the 6 and the time the computer receives this com- 
mand (d ). Some mechanization of the first order filter on the aircraft 
^ 1 

computer will generate one delay interval (T = T ) . See Fig. IV-11. 

Ct 

is the output of the digital low pass filter. 



One of the delay intervals, T 2 , can be eliminated if an analog low 
pass filter is implemented. However, in this case, we are losing the 
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option of changing the time constant of the low pass filter for different 
flight conditions in the digital computer. We assume that the low pass 
filter [as in DO-1 and SU-l] is mechanized on the digital computer and is 
a part of the discrete control system. The control configuration des- 
cribed in Fig. IV-10 is closed by an additional feedback loop--the pilot, 
who is sensing 6 (the pitch angle). The pilot T s transfer function in- 
cludes a pure delay (At £0,3 sec), and together with the delay caused 
by the digital system, this loop may be unstable. A current research 
done by Stapleford [STA-lJ shows that by lowering the sampling rate to 

co = 10 to 15 cps, the pilot rejects this sampling rate before instability 
s 

occurs. The reason for this rejection is the roughness of control. 

The simulation results for different flight conditions and differ- 
ent sampling rates are analyzed and evaluated as follows: 

(a) Flight condition No. 1 (zero altitude, Mach 0.19). The be- 
havior of the aircraft is plotted in Fig. IV-12. The first 
plot is the response of the free (uncontrolled) aircraft to 
an elevator step-input. As shown, the response is slow and 
unsatisfactory ( TRP £ 1.0). Fig. IV-12 shows that after the 
closed loop pole relocation, the response to 6^ is fast but 
with an excessive pitch overshoot (TRP £ 0.27). By filtering 
the pilot stick input, a good time response was obtained 

(TRP £ 0 . 1 ) . 

The time response for lower sampling rates is plotted in Fig. 
IV-13 and Fig. IV- l4. As far as the TRP criterion is considered, 
the limit on the sampling interval is in the vicinity of T = 0.1 
(Fig. IV- 13 ). As seen in Fig. IV-l4, the slow response (TRP 
> 0 . 25 ) is caused by the two delay intervals, and T^. 

( b) Flight conditions No. 8 (zero altitude, Mach 1.2). The free 
flight behavior of the F-H aircraft, in this flight condition. 
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is characterized by a fast, short period oscillation with a low 
damping. The time response of the controlled aircraft to a 
filtered stick input for different sampling rates is described 
in Figs. IV- 15, IV-l6, and IV-17 . 

The time response parameter is within its specifications, but 
the TRP is not a suitable criterion for this discrete case. 

As seen in Fig. IV-17, the time response of the pitch rate be- 
tween the sampling points is essentially a free oscillation. 
Therefore we may conclude that, for this flight condition too, 
the acceptable limit on the sampling rate is in the vicinity 
of T =s 0. 1 sec. 


E. SUMMARY 

1. To preserve acceptable time response, there is a minimum value 
of the sampling rate due to the delays introduced by the sampling of the 
pilot 1 s filtered input. This minimum sampling rate depends on the required 
speed of response and on the time constant of the pilot 1 s stick filter. 

2. A sampling rate selection study was made for controlling the 
short period motion of the F-H airplane. The maximal sampling interval, 

T, was found to be in the vicinity of T = 0.1 sec; the rigid body 
short period was 10.5 sec and 0.476 sec for Flight Condition No. 1 and 
No. 8 respectively, and the first bending mode period was 0.2 sec. 
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ss CONTROLLED AIRCRAFT RESPONSE TO FILTERED STICK INPUT AT t - O 
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FIG. IV-12 


THE F-H SHORT PERIOD TIME RESPONSE TO A STEP 
INPUT. Flight Condition No* 1. 



FIG. IV-13 


THE TIME RESPONSE OF THE F-H SHORT PERIOD TO 
A FILTERED STICK INPUT. Flight Condition No 
1. T = 0.1 sec; TRP = 0.18. 



FIG. IV-14 THE TIME RESPONSE OF THE F-H SHORT PERIOD 

TO A FILTERED STICK INPUT. Flight Condition 
No. 1. T = 0.2 sec; TRP = 0.3. 
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V* SAMPLING AND ROOT MEAN SQUARE RESPONSE TO NOISE 


One of the principal objectives of the aircraft autopilot is to 
reduce the influence of random gusts of wind. The problem of gust allev- 
iation is especially severe during high speed flights at low altitudes. 

If a discrete control is used, there is a finite time interval between 
corrective commands. During this time interval, the aircraft’s average 
random motion due to external noise is increasing. If the time interval 
is too long, i.e., the sampling rate is too low, the rms value of some 
of the aircraft states may exceed acceptable boundaries. This places 
an upper limit on the sampling interval. 

The first systematic approach to sampling rate selection as it 
relates to rms response was devised by Berman [BE-l] in 1973. Berman 
discusses measurements contaminated *by white noise at the sampling 
instant. Then he propagates the state covariance matrix of an aircraft 
disturbed by white noise until one of the covariants goes beyond a pre- 
determined limit. The elapsed time is the sampling interval. Using this 
interval, Berman calculates the optimal (Kalman) observer. Controls are 
predetermined by a pole placement. 

However, the purpose of an optimal observer is not only to recon- 
struct the unmeasured states. In addition, it makes a better estimate 
of the states than could be obtained from noisy measurements alone. 

The uncertainty of the measurement taken at the beginning of the sampling 
interval is only known after the optimal observer has been determined. 

But the optimal observer can be calculated only if the sampling interval 
is given. For this reason, Berman’s approach would yield a higher value 
of the covariance than would result from using an optimal observer. 

We will show a systematic way to select a maximal sampling interval 
related to the rms response. 
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A. PRELIMINARIES 


The behavior of a continuous system driven by white or colored 
noise is well understood and is described in various sources [BR-1, 

BR-4, KW-l], A closed loop system, with a perfect observer and a fast 
controller, may reduce the influence of an external colored noise to 
any degree of acceptable rms response provided that a proper (quadratic) 
criteria is assumed and that no restriction is imposed on time response. 

In the case of nonminimum phase systems, or systems with a limited 
amount of control available, the rms response will depend on dynamic 
characteristics of the closed loop system and on correlation times of 
the colored noise. Discretely controlled continuous systems, driven 
by a continuous noise, white or colored, are more complicated. It 
will be proven that if an unlimited control were available, there is a 
well defined limit on the noise alleviation. This limit depends on the 
sampling time. 

The objectives of this section are as follows: (l) to outline 

the basic relationships; ( 2 ) to demonstrate the limitation on noise 
reduction due to sampling. 


A-l Basic Relationships 


A continuous plant, F, driven by a colored noise w, may be 
formulated as 



(5.1) 


where T) h> n(0, q), white noise, and w is n order Gauss-Markov 

tJ 

process for the case of a discrete control; the equations are: 


r * 

X 

li 

i — i 

e 

X 

- w - 

i+1 ! 

w 


i 


+ [r x ] u. + [r 2 3 T| 1 , 


(5.2) 
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where -> N(o, Q d ) • 

As described in Chapter III, for practical reasons, we shall choose 


r 


2 


= I 



and consequently, Q_(t) will be given as 

a 


Q d (T) - <I(T)X(0)$(-T) T 


T 

+ J *(T - t)G 2 QC2<6 T (T - T)dT . 

O 


( 5 . 3 ) 


Assuming X(o) =0, = Cx^, the rms response (or the steady state 

of the GaussrMarkov sequence) is the solution of Bryson and Ho fBR-l], 

X - (<r + r^xc® + r 1 c) T + Q d (5. if) 


where X is the covariance matrix of the states; the square roots of 
the main diagonal are the rms responses of the states. $ + T^C is the 
closed loop transition matrix, but not necessarily the optimal one. 

In this case the state of the wind w is not fed back. Equation 5*^ 
assumes a perfect knowledge of the states of the plant. 

If the measurements are highly contaminated by noise, the rms 
response of the system can be calculated by augmenting the plant equa- 
tions, (5.l), by observer equations and solving for X from (5*4) for 
the augmented system. 

In Section V-2, we will use the optimal control and optimal filter- 
ing approach which was described in detail in Chapter III. 

A -2 Basic Limitation on Noise Alleviation Due to Sampling 

Contrary to the continuous control systems, the discrete commands 
are executed at discrete intervals. During these intervals, the system 
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is exposed to disturbances. The best way to reduce the noise influence 
is to try to use inputs which will bring the system to zero as quickly 
as possible. But there is a limit to the noise alleviation no matter 
what control we use. This will now be shown. First, some general 
results will be derived. 

(a) A linear discrete system controlled by a state variable feed- 
back, can be reduced to a zero state at most in n steps 
(n- order of the system). A complete proof, using the Caley- 
Hamilton theorem, is given in Kwakernaak [KW-l]. 


If 4 

c 


is the closed loop transition matrix, and 


x i+ i 


4 x* 
c 1 


(5.5) 


then the resuits are that (l) $ is a nilpotent matrix 

c 

(4 = 0), and (2) all eigenvalues of 4 are located at 

c c 

the origin. In classical control literature, this system is 
named a "deadbeat” system, or a 11 finite settling 11 time design. 


(b) Claim: Given a continuous system with zeros at the origin; 

controlled by an optimal discrete controller, then if all 
A^/b co (unlimited amount of control available), the system 
is equivalent to a deadbeat system, where A^ is the weight 
on the state x^; B is the weight on the control, assuming 
a single input system. 


Proof ; Using root sqiiare locus, the system (5-5) can be 
formulated as [KW-l] 

B + y(z" 1 ) T Ay(z) = 0 (5.6) 

where y(z) is the vector of the transfer functions, x(z) = 

y(z)u(z), (u scalar); and A is the weighting matrix of 

the states in the quadratic criteria. The transfer functions 

y.(z)* where x.(z) - y.(z) u(z) have the following structure 
J J J a 
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raj 

where z are the zeros corresponding to the state j, 
and P (z) is n-order characteristic polynomial of the 
open loop system. For physical systems, £ n. Using 

(5.7), Eq. (5.6) can be formulated as 



(5.8) 


m j 

or, by multiplying the numerators and denominators by z 
and rearranging 

J [BP n (z mJ ’ 1 ) p n (z) + = 0 

as A^/b <», the poles of (5.7), which are the poles of the 
closed loop system, approach the zeros of the open loop system 
(located at the origin), and therefore become a deadbeat sys- 
tem. 

Using these results, we can formulate the following theorem. 

Theorem : For a discretely controlled continuous system, the rms response 

to white noise, for any control, is equal to or greater than: 

T 

X = + ... + « c Vc + Q d • (5.10) 

If all the closed loop zeros are located at the origin, the equal- 
ity holds. 

Proof ; The closed loop system, (5»5)> driven by white noise, is 
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<5 x. + rw, . 

ci i 


(5.ii) 


X 


i+1 


From (5.1l) we obtain 

x = <r n x. + ® n *rw + . + i_a w (5*12) 

i+n c i c i 1 

T 

but $ n = 0. Multiplying (5.12) by x._. and averaging yields 
c i+n 

T 

E( Xi X^) = X a $ n-1 Q d $n_1 + + Q d • (5.13) 

i>n 

Equation (5.13) gives the lower limit on X. Note: (l) <$ c = * C ( T ); 

(2) X = X(T,Q); (3) If not all the closed loop zeros are located 

at the origin, or A /b is limited, then / 0, in this case, 

$ n X x T $ nT > 0 and the rms response to white noise is greater than 
i i 

(5.13). 

Example: A first order system, driven by a continuous white noise and 

controlled by a discrete controller, e.g., the following system 


x = ax + gu + w 


w N(0, q) 


aT 

x. _ =* e x. + Tu + w 

i+1 1 i 1 


w 


± -* N(0, q d ) (5.14) 


q = \ e 'qe 'dT 
d 


at aT = S. (e 2aT -1) 


2a 


Closing the loop of (5.l4) by Tu = bx^ yields 


X i + 1 


aT 

( e + bj x^ + w 1 


N[0, — (e 2aT - l)]. (5.15) 


The rms response of the closed loop system to white noise is the 
steady state solution of the Gauss^Markov sequence (5*13)* 
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o 

Definition: at. = q - 

x i Xl 

by the solution of ( 5 «^+) 


The steady state value of q Xi is given 


q , 2aT 

(e 


Sc, 


2 a 


1 ) 


(5.16) 


ss 


/ aT w \ 2 

(e + b) 


Note: for T = 0 (continuous control), q = q/-2(a+b), a + b 

s s 

< 0* Some properties of (5.l6) are: (a) for b = 0 or T -* 

q x , -q/ 2 a which is the open loop steady state covariance of x. 

(b) for a given T, the deadbeat b is 


aT 

e 
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Equation ( 5 . 17 ) shows that even for any control available, q 

^min 

is always finite, and zero response to external noise can be 
achieved by limiting T -> 0 only, T 0 =£> q x ^ -> 0. 

This analysis is strictly valid on sampling points only. If 
the open loop is stable and the output of the controller is a deadbeat 
system, then this analysis could be extended to intersample points also. 


B. DETERMINATION OF THE MAXIMUM ALLOWABLE SAMPLING INTERVAL 

We will analyze in detail the influence of various factors which 
determine the rms response of a discretely controlled continuous system. 
Using the expression for rms response calculation derived in Chapter III, 
we will consider the influence of 

1. Different correlation times of external disturbance; 

2 . Optimal control vs pole placement; 

3. Varying accuracy of measuring instrument. 
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A clearer understanding of these factors will help us to select the 
proper sampling rate. Our calculation will be carried out for the F-H 
aircraft example, flying in the most gusty environment. 

After detailed analysis in Sections V-B-l, V-B-2, and V-B-3, we 
will summarize the method of selection of the sampling rate in Section 

V-B-4. 


B-l Relation Between Correlation Time of The Colored Noise Disturb - 
ance and the Sampling Time 

For the rms calculations we will use a simplified model of the 
short period dynamics: 
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1], T) -» N(0, l). 


The wind gusts are modeled as a colored noise disturbance. "Colored" 
means that the power at the wind gusts is not equally distributed over 
all frequencies. If, as in our case, a first order Gauss-Markov 
process is used for the stochastic wind gust model, then the correla- 
tion time determines the spectrum of the wind gust. More power is 

concentrated in the lower frequencies than in the higher frequencies. 
The rms value of the wind is related to the integral of the power 
spectrum and is designated by a . The mathematical process which 
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generates this colored noise is based on a first order filter, with the 

time constant t , which filters a white noise with a power spectral 
2 / W 

density 2a /r - 
' w 

Modeling the disturbance as a colored noise instead of as a white 
noise has three advantages: 

(a) It is a more realistic model from an engineering point of 
view ; 

( b) The white noise has a weak mathematical definition [KW-l]; 
moreover, its total power is infinitely great* 

(c) In the short period example, one of the measurements is the 
acceleration. The acceleration is related to OL ^ which is 
combined from the angle-of-attack and wind gusts. If a white 
noise is used as a model of the wind disturbance, then the 
measurements are contaminated, not only by the measurement 
noise, but also by the process noise. This makes the calcu- 
lation of the optimal observer much more difficult [BR-l] # 

The colored noise can be characterized by its spectral distribution, or 
by its correlation function. The correlation function C^(T f ) indicates 
how the disturbance at time t relates to a disturbance at time 
t -h t 1 . The various relationships can be visualized in Figs. V-l and V-2. 



FIG. V-l THE SPECTRUM OF Iw I . 

1 g' 
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FIG. V-2 THE CORRELATION FUNCTION OF w . 

g 


To calculate the rms response to a vertical gust, we will compare 

the influence of different sampling rates and different correlation 

times. The most meaningful comparison can be made if the gusts have the 

2 2 

same intensity , i.e., a - constant. In this case, a , which is the 
variance of the wind gust, is related to the intensity of the wind. The 
energy distribution is schematically described in Fig. V-3 ( log scale). 



FIG. V-3 ENERGY DISTRIBUTION OF WIND WITH THE SAME 
INTENSITY AND DIFFERENT CORRELATION TIMES. 
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The area under each curve is related to the total kinetic energy of the 
wind and is constant. 

Using the methods from Chapter III, the rms response of q and 0^ 
were calculated for different sampling times and different correlation 
times. 

The numerical example we are investigating in this Chapter is the 

flight condition No. 8. of the F-H aircraft described by Borow [BO-l] and 

Sutton [SU-1]. This is a flight in zero altitude and at Mach 1.2. It 

is usually assumed [cf ET-1 ] that the wind disturbance is “frozen" 

in space and the velocity of flight determines the correlation time ? . 

w 

Furthermore, depends on various other factors such as the terrain 

over which the aircraft is flown (zero altitude). Borow and Sutton 

each give the average values of a and t as o' = 12 ft/sec and 

w 

T = 0-5 sec. During our investigation, we will vary r only. The 
w w 

total rms response is directly related to a and therefore it will not 
be necessary to analyze the influence of various o' s on Q!^ and q. 

The rms response depends also on the weighting matrices in the quadratic 
cost function. In Chapter IV we described how the weights on the states 
are chosen. It was shown that the selection of the weights represents 
a compromise among such factors as rms response, time response to a 
pilot input, bandwidth of the actuators, and of the closed loop system. 
Taking into account these considerations, the selected weights and the 
corresponding closed loop roots are summarized in Table V-l: 


Table V-l 

LOCATION OF THE CLOSED LOOP ROOTS, FLIGHT CONDITION No. 8 


Closec3^\^^Sampling 
Loop Roots 

T = 0.005 sec 

T =0.05 sec 


z-plane 

0. 929±j 0. 045 

0.41 ± j0.22 

0. 065±j0. 15 

s-plane 

—14 . 5±j9. 6 

-16 ± 12.0 

-17±jl2.5 


Open Loop Roots: s =-2,5 ± 13.5 

Weights: A /B = 0.12 (sec S A = 11.0 [01. 

q °t/b 
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The results of the rms response calculations are plotted in Figs, 
V-4 and V-5 as a function of various sampling times and various correla- 
tion times, q is the rms response of the pitch rate. Reduction of the 
pitch rate of the F-H example when the aircraft is flying at zero 
altitude, is important for precise aiming. 0^ is the rms value of the 
total angle of attack, and includes the angle of attack generated by 
the vertical gust, a is directly related to the lift and therefore 
to the vertical acceleration, which in turn determines the pilot's 
comfort. Figure V-4 demonstrates how the rms response increases as 
the sampling interval increases. Figure V-5 demonstrates that a system 
disturbed by colored noise is better able to alleviate noise with a 
longer correlation time and identical intensity. 

Certainly the designer cannot choose the correlation time of the 
vertical noise; this is an empirical, measured quantity. But these 
figures show that it is important to identify the correlation time of 
the disturbance and not only the rms value of the gust. For the same 
intensity of the vertical gust, the designer can choose a longer samp- 
ling interval if he knows the correlation time of the external noise. 

B-2 Reduction of RMS Response by Optimal Control Compared to the 

Classical Approach . 

As explained in the previous section, there is a definite rela- 
tionship between sampling time, correlation time, and rms response. 

Gust alleviation can be achieved by classical control methods, essen- 
t ially by pole relocation, but maximal reduction can be obtained by 
using existing knowledge of noise characteristics. Colored noise is 
characterized by its correlation time constant. For longer correla- 
tion times, we may make a more accurate prediction of future behavior 
of the disturbance and consequently we may apply better feedback control 
in order to reduce the disturbance influence. We will describe in more 
detail how the optimal controller uses the information on behavior of 
w^ in the following equation, repeated from Eq. (5.2): 
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S i+1 


+ r 1 u ± + [r 2 ] ti. ^ -» n(o, q ) 


The matrices <f. T, } P and the covariance matrix Q. have the 
7 1 2 a 


following form: 


*1 *12 


Q<l ll 


Si V 


r 3 " 1 


Using the separation theorem, the optimal controller 


u. = [C C 
x x 


■’ [:]. 


is calculated from the deterministic part of (5.2). Furthermore, it 

could be proven that C is unaffected by the extended state of the 

x 

external disturbance w [HAL-1 ], i.e., C is not a function of 

g x 

and $ . In other words, C feedback gains of the system states 

i ^ x 

are unaffected by the state of the disturbance, but not vice versa. 

The feedback gain of the external disturbance depends on the 
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structure of the system. This indicates the main advantage of the 
optimal design. Even if the closed loop roots are relocated by other 
methods to the same position, the feedback gain always reduces the 

rms response. 

The rms response for the augmented short period mode, including 

wind as a feedback state, was calculated for different sampling times 

and different correlation times and is summarized in Fig. V-6. Figure 

V-6 describes the rms response at the pitch rate to a vertical gust 

for two different designs: (a) C / 0 is the optimal control approach; 

w 

( b ) C w =0 is equivalent to a pole placement design, the state of the 
disturbance is not fed back. Figure V-6 clearly shows that the optimal 
approach will be preferred, especially at low sampling rates and for 
shorter correlation times of the external disturbance. 

B- 3 RMS Response for Imprecise Measurements 

In this section we will analyze the influence of measurement 
noise on the rms response to the vertical gust with reference to the 
F-h example. We have no precise knowledge of the states. Only q 
(pitch rate) is measured directly. The other two states, OC^, (angle 
of attack) and (wind gust) are estimated by an observer. For our 

illustrative example, we will not use a reduced order observer even 
if q is directly measured. 

The measurements are made by an accelerometer (n ), and a 

a 

rate gyro, q. 
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The power spectral density matrix R of the additive white noise 
is calculated by the following method. There are two major sources of 
measurement noise. 

1 . The A/D converter resolution inaccuracy . This error is 

regarded as a white noise independent of the sampling rate. 
Borow [BO-1 ] assumes the power spectral density of this noise 
to be 0.1$ of the full range of the measuring instrument. 


2 . 


The Instrument noise . Borow considers this noise to be colored, 
with a correlation time of about r = 0.01 sec and an rms 
level of 0 . 25 $ of the full range of the measuring instrument. 
This continuous colored noise can be approximated as a dis- 
crete white noise by a procedure suggested by Bucy and Joseph 
[bU- 1 ], where the discrete covariance of the sampled noise is 


given by 


R „ 2 1 + *' TA 


The full range of the instruments is approximated in Borow, 

for the accelerometer n = 10g, and for the rate gyro, 

max 7 

q = 1 rad/sec* 
max 

The total measurement noise is derived from these two sources. 
The equivalent power spectral density matrix of the discrete noise is 
given in Table V-2, 


Table V-2 


THE POWER SPECTRAL DENSITY OF THE MEASUREMENT NOISE 


T 

R 

n 

R 

q 

0.005 sec 

2 , 

2.6 ft /sec 

-5 2/2 

2,6 X 10 rad /sec 


2 , 4 

-6 2 7 2 

0.05 sec 

0.8 ft /sec 

8.0 X 10 rad /sec 



-6 2> 2 

0. 1 sec 

0.8 ft /sec 

8.0 x 10 rad /sec 
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The rms response of the closed loop system was derived in 
Chapter 3> Eq. (3*86). We repeat it here as a solution X of 

x - m = O + rc)(x - p)(4 + rc) T . 

Obviously for larger R, the rms response of the closed loop systen 
will be larger. 

In Fig. V-7 the rms response of the pitch rate and the total 
angle of attack are plotted for different levels of the power spectral 
densities of the measurement white noise. The price of the measuring 
instrument is related to its accuracy. Precalculated figures, such 
as Fig. V-7, can help the designer trade off between the sampling rate, 
the external noise reduction, and the cost of the measuring instruments. 

B-4 Selection of Sampling Rate for the F-H Example 

The properties of the external disturbance were given in Section 
V-B-l . They are a = 12 ft/sec, =0.5 sec. In Borow and Sutton, 

no limits were given on the acceptable rms response. The general 
requirement is to reduce the gust influence as much as possible. 

Using the optimal control and the optimal filtering approaches, 

the rms responses were calculated and plotted in Fig. V-8. The dynamic 

properties of the closed loop are achieved by proper weightings in the 

quadratic cost function (Ch. IV), keeping in mind the bandwidth of 

the actuator and the time response to a step input. From Fig. V-8 we 

can see that for sampling rates higher than co =20 cps (T = 0.05 sec), 

s 

little improvement in gust alleviation is realized. Therefore, when 
gust response is considered, the sampling rate of the F-H short period 
mode should be in the vicinity of = 20 cps. 
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C. SUMMARY 


1. For a discrete controller, the lower limit on the rms response 
to a random disturbance is a function of the sampling frequency, the 
rms disturbance level, and the correlation time of the disturbance. 

2. By using a simplified model, and a stochastic model of the 
external disturbances, the designer may determine the lowest sampling 
rate which will yield an acceptable rms response, 

3. Furthermore, by plotting graphs such as Fig. V-7, the designer 
may choose measurement instruments, which are economical with respect 

to rms response. 

4. We conclude that for our principal example, F-H flight condi- 
tion No. 8, the sampling rate should be between 10 and 20 cps, in 
order to reduce effectively the rms response to gust. Any increase in 
the sampling rate past 20 cps results in a very small improvement. 

5. An important consideration is that we are using linear models 
and therefore the weights on the states should be chosen, keeping in 
mind the bandwidth of the actuator. 
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VI. SAMPLING TIME AND SENSITIVITY 


In the linear discrete system theory, if a perfect linear model 
is given, there is no lower limit on the sampling rate, i.e., one can 
obtain closed loop roots in an arbitrary location in the z-plane for 

any (o (sampling rate). A completely different problem is the sensi- 

s 

tivity of the closed loop poles to variations of different parameters 
of the assumed model. For example, in the F-H short period extended 
model with the nearly undamped bending mode, the measured signals of 
the accelerometer and of the rate gyro are contaminated by the bending 
oscillation. Those combined measurements are fed to the observer 
which reconstructs the desired states of the system. As shown later, 
the optimal compensator generates a notch filter which filters the 
unwanted bending frequency. But if the actual bending frequency differs 
from the assumed model, the incoming bending frequency misses the notch 
and is fed back to the elevator as a positive feedback. Potentially, 
this mechanism can destroy the stability determined using a perfect model 
of the system. The purpose of this chapter is to analyze this situation 
and to find a way to reduce the sensitivity to bending frequency variation. 
We will investigate the following aspects of this problem. 

In Section VI-A we will show that in a linear control system, which 
includes a nearly undamped and uncontrolled frequency, the optimal com- 
pensator generates a notch filter. The depth of the notch is directly 
proportional to the amount of the sensor pickup of the unwanted frequency. 

In Section VI-B we will investigate the sensitivity of the short 
period example to the bending frequency variation. We will show that the 
sensitivity increases for lower sampling rates. 

In Section VI-C a method which desensitizes the system will be 
shown. The method is based on relocation of the observer error poles 
and of the closed loop bending poles. The primary goal is to increase 
the width of the notch filter rather than to control the bending mode. 
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In Section VI-D we will show how the rms response was changed alter 
the desensitization of the optimal system (as described in VI-C), It 
will be shown that in this case, the rms response will increase slightly. 
This is the penalty we have to pay for a less sensitive system. 

In Section VI-E we will show that there is a definite relation be- 
tween the location of the compensator poles and the sensitivity of the 
closed loop system. A system with a more stable compensator is less 
sensitive to variations of parameters. A theoretical proof will be 
given. 

In Section VI-F we will look into the possibility of filtering the 
the bending frequency by sampling at the same rate as the bending oscilla- 
tion. It will be shown that this approach may cause a hidden excitation 
of the bending mode. This instability cannot be detected by a regular 
analysis in the z-plane. 


A. FORMULATION OF THE PROBLEM 


A- 1 Basic Configuration 

The basic configuration being investigated is described in Fig. 

VI-1. 



FIG. VI-1 THE BASIC SYSTEM FOR SENSITIVITY INVESTIGATION 
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The oscillating modes, in our case the structural dynamics, are 
part of the system. The primary purpose of the compensator is to control 
the rigid body dynamics, but it must be done in such a way that the bend- 
ing modes are adequately damped. In the basic system (Fig, VI-1), we 
assumed that there is no direct coupling between the plant and the oscillat- 
ing modes. This is the usual approximation found in several sources [BO-1], 
[BL-1], [SU~1], for modeling the rigid body dynamics and the structural 
modes* There will certainly be an important coupling if the closed loop 
bandwidth of the plant will be too near to the bending mode frequency. 
However, usually the bending mode frequency is higher than the closed 
loop frequencies* In our short period example, the first bending mode 
(the slowest), is located at = 25 rad/sec, while the fastest short 

period open loop frequency is 13,5 rad/sec. The closed loop short 
period mode frequency is about 10 rad/sec (flight condition No* 8, see 
Ch* IV), 

Using the results from Chapter III, the different blocks of Fig, 

VI— 1 will be described in more detail in Fig* VI-2, where the following 
definitions have been made: 











As described in Chapter III, the compensator transfer function is based 
on 


= + K(v - Hx )\ 


• observer 


h + i * + ru i 


therefore , 


u. 

1 


u. - Cx. 

1 l 


C[(I-KH)(4>+rC)-Iz] 1 Kzy. 


controller 


compensator* 


In this simplified system, the oscillating bending modes are separated 
from the system, and are only slightly influenced by external disturbances. 
Sutton [SU-1] and Blakeiock [BL-1] assume that the bending modes are not 
directly influenced by the external disturbance* Their assumption is 
that the bending mode is excited by the elevator input only* Borow [BO-1] 
uses a more elaborate model and assumes some direct influence of the 
external disturbance on the bending mode. We are using the approach of 
Borow, which is more realistic* 
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A- 2 Optimal Compensator as a Notch Filter 


The compensator transfer function (6.3) can be more explicitly 
written as 

T t 

C adj [ (I-KH) ($+rC)-Iz]Kz 

u(z) = • (6.4) 

|[(I-KH)($+rC)-Iz] | 


Using a well known theorem from linear systems theory (e.g. , [CH-1]), 
the zeros of a linear system are unaffected by a state variable feed- 
back and therefore (6.4) can be simplified to: 


C adj T [<t - Iz ]Kz 

u(z) = y(z) . (6.5) 

| [(l-KH)(<D+rc) - Iz]| 


Note that in the case of the compensator, the state variable feedback is 

derived from gains of the observer and of the controller. Furthermore, 

using the model of (6.1), where and 0^ aI * e uncoupled, and assuming 

that no attempt is made to control the bending modes, i.e. # C =0. 

b ' 

the numerator of Eq. (6.5) obtains the form; 



"(« -Iz )" 1 

0 



K 

a 


a 


|$> -Izj U-Izl 




0 


a ' 1 b 1 

j 


A. 


( 6 . 6 ) 


and the simplified compensator transfer function will be 


where adjf ] = cofactor matrix of [ ]. 
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u(z) = 


$ -Iz 
b 


C adj T [® -Iz] K 
a a az 


— ,y(z) 


ri-K h 
a a 


K H, 
a b 


t / 1 XT T— TC H 

— i'l, ii d-iv n 

b a b b 


rtf +r c 

a a a 

cT 


~Iz 

0 ^ 

- r* p 

b a 

flv J 

1 

- 0 

Iz- 


(6.7) 


Johnson [JO-1], using the relation of (6.6), points out that the zeros 
of the compensator (represented by - Iz|), consist of a notch filter 

which is located exactly at the bending frequency. But a proper notch 
filter is a combination of second order zeros and second order poles 
[CA-2], as 


notch filter 


(« - r 1 . J “‘ r )(* - r ie ' J “' T ) 

(. - y J “* T )(z - r 2 e-3“^ 


( 6 , 8 ) 


Here, is as near as possible to (co^, the tuned frequency). 

The sharpness of the notch is related to the ratio r /r^. r i is relatvCi 
to the damping of the open loop bending mode and cannot be changed. By 
manipulating the controller gains C and the observer gains K, the 
notch filter poles can be shifted to a suitable position. However, 
Johnson, who uses a pole placement approach, considers a completely 
different criterion when observer error gains are chosen. Their cri- 
terion is a proper time response of the whole closed loop system. 

Usually the time response of a controller— observer design is unaffected 
by the observer error g^ins. However, Johnson doesn’t feed the input 
commands to the observer; therefore the system response is significantly 
influenced by observer gains. As indicated by Eq. (6.2), the pilot input 
and the measurements are fed to the compensator. Except for sensitivity 
problems, the observer error gains are chosen by the optimal design from 
the standpoint of noise statistics. To obtain a better understanding 
of the filtering action of the compensator, we will investigate the 
behavior of the compensator for the extended model of the short period 
example. 
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As pointed out at the end of Section VI-A-1, several references 
assume no direct influence at the external disturbance on the bending 
mode. From a sensitivity point of view, that is a highly undesirable 
model. The optimal design procedure in this case considers the bending 
mode to be an undisturbable mode and generates zero gains for the ob- 
server's error. Consequently, the damping at the observer's error of 
the bending mode is nearly zero. We are using the approach of Borow 
[BO-1], which assumes a direct influence of the external disturbance on 
the bending mode. Hence, the optimal procedure generates a higher 
damping of the observer error of the bending mode. See Table VI— 1. 

Table VI-1 

NOTCH FILTER POLES-ZEROS LOCATION. NOMINAL CONFIGUR- 
ATION, FLIGHT CONDITION No. 8 


Open loop bending mode: 
(notch filters, zeros), 
T = 0.05 sec) 

z 

= 0.33 

± jO.94 

s-plane equivalent : 

s 

= -0.5 

± J25 

Observer error bending 
mode poles, T = 0.05s: 

z 

= 0.32 

± jO.91 

s— plane equivalent: 

s 

= -0.98 ± J24.6 

Compensator bending mode 
poles, T = 0.05s: 

(notch filter poles) 

z 

= 0.25 

± j 0. 91 

s-plane equivalent: 

s 

= -1.3 

± J26.5 


The filtering properties of the notch are characterized by the ratio 
of the open loop bending poles which are the notch zeros, and the com- 
pensator bending mode poles. As we can see from Fig. VI-3, the notch is 
generated by superimposing the frequency response of the poles and the 
zeros of the notch filter. 
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FIG. VI-3 NOTCH FILTER GENERATED BY THE OPTIMAL 
COMPENSATOR. 


The amount of the bending mode pickup, depends on the loca- 

tion of the measuring instruments with respect to the bending mode. If 
the accelerometer and the rate gyro are located exactly on their neutral 
points as per Fig. VI-4, then there is no direct measurement of the 
bending mode (H^ = 0), and the numerator bending zeros (in Eq. 6.7) 
are cancelled by the denominator bending poles. 


h = 0 h = 0 
r a 



FIG. VI-4 THE NEUTRAL POSITION OF THE ACCELEROMETER (h ) AND OF THE 
RATE GYRO (h r ) WITH RESPECT TO THE BENDING MODE. 
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It is highly impractical to assume that the measuring instruments can be 
allocated to the neutral points. The neutral points (H^ = 0) are 
changing during the flight. Moreover , the neutral points for higher order 
bending modes (not considered here) are not in the same positions. There- 
fore, we will analyze the behavior of the notch as a function of H^. 

In Fig. VI-5 the compensator bending pole location is traced as a func- 
tion of the accelerometer position, assuming that the rate gyro is in 
a neutral position (h^ =0). As can be seen, the frequency increases, 
but the damping of the poles is nearly unchanged. Looking at Fig. VI-3, 
as the frequency increases, the frequency pattern of the poles moves to 
the right hand side and the depth of the notch increases. 

It should be emphasized that in Fig, VI-3, only the notch response 
is plotted. The total frequency response of the optimal Compensator in- 
cludes all other factors of (6.7), 

B . SENSITIVITY OF THE SHORT PER I 0D MODE 


Based on the previous section, which shows how the optimal com- 
pensator filters the unwanted frequency, we will investigate the sensi- 
tivity of the system to a variation of this frequency. 

B— 1 Sensitivity Definition 

The sensitivity of the closed loop system is defined as a change 
in the location of the poles due to the parameter variations of the 
system. The closed loop system includes the plant: 


x. , ~ $ x. + T u. + w. 

l+l pi pi i 


y. = Hx. + v. , 
l l l 9 


(6.9) 


and the compensator: (repeated from Eq. 6,2) 
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FIG. VI-5 


RELOCATION OF THE COMPEL 
FUNCTION OF THE ACCELEROfl 


OPEN LOOP BENDING POLE 



h 

r 


O, a neutral position 
of the rate gyro 


h^ = nominal position of the 
accelerometer 




x. = x. + K(y. - Hx. ) 
1 1 1 l 


x 


i + 1 


u, 

l 


$x_ + ru. 

i i 

/N 

Cx, . 

l 


If <r> = $ and r = P , the assumed model in the observer, is completely 
P P 

identical to the plant, then we will define the observer T s error x as 

~ A 

X = x - x # 

Combining the compensator equations with the plant equations of 
state, we get 


x = ($ + rC)(x. - KHx. - Kv ) 

l+l l li 

x = 0(1 - KH)x. + OKv. - w. . 

l + l l li 

The characteristic equation of the system (6.10) is: 

U+rc-iz -($+rc)KH ( 


= o. 


0 $ (I-KH)-Iz 


( 6 . 10 ) 


( 6 . 11 ) 


The closed loop poles, for a perfect knowledge of the plant parameters, 
are the poles of the closed loop plant and the observer error system: 

|3)+rc-lz [ Jo (I— KH)—Iz j = 0. (6,12) 

The zeros of expression (6.12) are the poles of the total closed loop 
system, (6.10). But by inspecting (6.2) and (6.10, we see that the 
time response of the closed loop is only influenced hy the zeros of 
|®+rc-lz| = 0. These are the poles of the transfer function u i y^ . 

It could be shown that the zeros of (6.12) are the poles of the transfer 
f unct ion v\ y^ , 
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If $ ± $ and T ?* T, then the plant-compensator system becomes: 

P P 


X = Ox. + r Cx, + w. 

i+1 pi pi i 


y. - Hx. + v. 

J i 11 


(6.13) 


x . = (I - KH) (0 + rC)x + Ky . 

1+1 i i+l 


The last equation can be rewriten as 


x. , = (l-KH)(0+rc)x. + KH(0 x. + r Cx. + w,)+ Kv. 

l+l ' 7 7 l pi Pi 1 1 


(6.14) 


The characteristic equation, for the case of an imperfect knowledge 
of the plant, will be 


0 -I 
P z 


r c 

p 


KHT (l-KH)( 0 +rc)+KHT C-I 

P P z 


= O . 


(6.15) 


We may still define x, the observer error, but separation of the poles 
no longer exists. The poles of the closed loop system are the eigen- 
values of (6.15). 

The sensitivity can be defined as the relocation of the eigenvalues 
of the system (6.15) with respect to the perfect system, (6.11). 

The pole's location on the z-plane does not provide us with a 
clear physical interpretation of the behavior of the system, particularly 
for poles located near point 1 on the real axis of the z-plane. Therefore, 
as an aid to evaluating the performance of the closed loop system, an 
equivalent damping, , based on sampling of a continuous system will 
be defined. See Fig. VI-6, 


o 
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z -PLANE 



FIG. VI-6 z- PLANE DESCRIPTION OF C*. 


z = 


z + iz 
a b 


= e 


e = idt 


■arc* 




T 1 


-l z b 
tan — 
z 

a 


( 6 . 16 ) 




-log 


(^ z a + Z b) 


{[ tan ‘ 1 ( r } ] + 


? * 


We may consider the equivalent damping ratio £ * as a damping ratio of a 
continuous sampled second order system* This definition is necessary 
since the discrete system, controlled by a zero order hold, behaves 
differently between the sampling points* 
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B-2 Sensitivity of the Extended Short Period Mode to a Bending 

Frequency Variation 

Because of the deep and sharp notch generated by the optimal com- 
pensator, we assumed that the system is highly sensitive to the bending 
frequency variation. We will investigate the sensitivity by varying the 
actual bending frequency by 10% from its nominal value, i.e., from 25 
rad/sec to 22*5 rad/sec. 

The model of the extended short period mode (flight condition No. 

8) was described in Chapter IV. Using the controller and observer design 
of Chapter IV, the sensitivity was calculated by solving (6*15) and 
(6.16). The behavior of the equivalent damping of the bending 

mode as a function of the sampling rate is plotted in Fig. VI-7. We 
can see how the sensitivity increases for lower sampling rates. 

The time response of the normalized bending mode, , to a wind 
gust impulse is plotted in Fig. VI-8 for = 0.02, 0.44. The 

simulation scheme is described in Appendix B. Note that, contrary to an 
input command, the impulse is applied only to the plant and not to the 
observer. Therefore, correcting commands are only activated after one 
sampling period. The main reason for the instability of the case 

= 0.44 (T =0.1 sec) is the undamped behavior of the observer 
error bending mode. In other words, the observer is unable to follow 
the plant. 

In the next section, we will describe a method for filtering the 
bending mode, even for an imperfect knowledge of the bending frequency. 

£. REDUCTION OF SENSITIVITY TO BENDING FREQUENCY 

VARIATION 

To reduce the sensitivity to the bending frequency variation, we 
have to increase the width of the notch filter which is generated by the 
compensator. As explained in Section VI-A, the width of the notch filter 
is related to the damping of the bending mode in the compensator. The 
damping may be increased by a modification in the control loop, or 
in both. This modification can be achieved by three basic methods: 
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FIG. VI-7 THE SHORT PERIOD MODE. Imperfect knowledge 
of the bending frequency. Q * vs w A* 
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FIG. VI-8 BENDING MODE TIME HISTORY. 10# error 
in V 
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(a) By a pole placement, 

(b) By using an automatic synthesis procedure, which will satisfy 
the quadratic criteria and will reduce the sensitivity. 

An appropriate computer program exists for continuous linear 
systems [HAD-1] but it is not yet available for discrete 
systems. 

(c) By a modification of the optimal discrete system. 

In our synthesis, we used the third approach. Stabilization of the 
perturbed closed loop system was achieved by weighting the bending mode 
state, x^. The gain of the observer was changed by assuming that more 
noise would act directly on the bending mode. This modification might 
have changed the dynamic properties of the system, but we will show that 
no significant change took place, 

C — 1 Desensitivity by Increasing the Damping of the Bending Mode 

We are invest igat ing the nominal design of the short period mode 
described in Chapter IV and in Section VI-B, In this design, no weight 
was given to the bending mode in the quadratic criteria. For desensi- 
tization of the system, we will damp the bending mode by weighting it 
The relocation of the poles as a function of the bending mode weight 

(A ) is given in Table VI-2. 

4 


Table VI-2 

CLOSED LOOP EIGENVALUES AS A FUNCTION OF WEIGHT ON 
THE BENDING MODE (T = 0.1 sec). 


A x 
_ X 4 . 

Short Period Mode 

Bending Mode 

0 

z = 0.137 ± j 0.280 

z = -0.781 ± j0.584 

-6 



10 

z = 0.14 ± JO, 275 

z = -0.669 ± j0.485 

-5 



10 

z = 0.16 ± j0.24 

z = -0.516 ± j0, 29 


Rigid body weights: A qy / B - °. 12 , A a y B = H. 

Flight Condition No, 8. 
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From the results In Table VI-2, we conclude that when A„ varies from 

-5 4 

zero to 10 , the dynamic properties of the rigid body are not signifi- 

cantly changed. The influence of this modification on the stability of 
the closed loop bending modes is described in Fig. VI-9, where §* has 

been plotted vs A x . As seen in Fig. VI-9, it is possible to stabilize 

4 

the bending mode for an assumed 10 percent error in the bending frequency. 
The time response of the bending mode to a wind gust impulse is shown in 
Fig. V-10. The gust impulse was translated via the initial value theorem 
to initial conditions formulation. The calculation of this time response 
was done by a digital simulation scheme described in Appendix B. . Note 
that the impulse is applied only to the plant and not to the observer; 
therefore, the first observer’s states and the first command appear only 
after one sampling interval. We can see that A > 10 stabilized the 

A 4 

system. Further stabilization, especially that of the observer error, will 
be made in the next section. 

C-2 Desensitivity by Increasing the Samiiing of the Observer Error Bend- 
ing Mode 

The observer error gains can be increased by introducing an arti- 
ficial external noise on the bending mode. This artificial noise is not 
included in the total rms response to external disturbances of the system. 
But the dynamic and filtering properties of the observer is no longer 
optimal in the sense that it minimizes the influence of all the noises 
acting on the system. This change in the total rms response for the 
modified system will be treated in detail in Section VI-D. 

By varying the factor G in the noise distribution matrix P , 

2 

the observer error gains were changed. The observer error poles corres- 
ponding to the bending mode vs the factor G 4 are given in Table VI-3. 

The equivalent damping of the closed loop bending mode (10% error on o ) 

b 

is plotted in Fig. VI-11. As seen for G 4 > 100, the bending modes are 
stabilized. The time history of the bending mode is traced in Fig, VI— 12. 
The primary stabilization effect is achieved by a damping of the observer’s 
error. When these phenomena are understood, it is possible to desensi- 
tize the perturbed short period mode. 
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FIG. VI- 12 BENDING MODE TIME HISTORY FOR VARIOUS NOISE 

MAGNITUDES. 10 °j 0 error in oj . co / go = 0.5, 

T = 0.1 sec. A Y/1 = 10 - 6 . b * b 
7 x 4 
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Table VI-3 


OBSERVER ERROR BENDING MODE EIGENVALUES AS A 

FUNCTION OF G, , THE ARTIFICIAL NOISE (T - 0.1 sec). 
4 


G 4 

z- plane 

Observer Error Poles 

s-plane Equivalent 

0 

-0.725 ± j0.572 

-1.0 ± J25 

100 

-0.6736 ± J0.508 

-1.5 ± J24.8 

200 

-0.622 ± jO.452 

-2.5 ± J24.5 


D. THE EVALUATION OF THE PERFORMANCE OF THE DESENSITIZED 
CLOSED LOOP SYSTEM 


The rm s response of the optimal system was evaluated in Chapter V. 

In this chapter, Section VI-C, we desensitize the optimal design by 
changing the weights in the cost function and by introducing an artificial 
noise on the bending mode. As a consequence of this modification, the 
control and the observer error gains were changed. The modified closed 
loop system is not optimal in the sense that it minimizes the rms 
response when the noise statistics and the original weight in the quadratic 
cost function are given. Yet, from a more general point of view, this 
modification would be a-, better engineering solution if the performance of 
the system didn’t greatly change. The performance criteria are mostly 
characterized by the closed pole location and the rms response. 

As has been shown in Table VI-2, the rigid body poles of the short 

period example were only slightly relocated. The modification we used 
— 6 

(A x =10 ) negligibly influenced the time response to a pilot input. 

4 

The filtered pilot commands (Ch. Ill) are fed simultaneously to the air- 
craft and to the observer. During the first second, the time response 
is dominated by the closed loop rigid body poles and not by the observer 
error poles of the bending mode. However, the observer error poles and 
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gains of the bending mode were changed (Table VI-3), and thus the steady 
state rms response properties were also changed. Therefore, it is 
necessary to find out how the modification of the closed loop system 
influences the rms response to external disturbances. This investigation 
will be conducted by considering the total rms response of the modified 
closed loop system to all noises acting on the system. 

The relation for the states covariance matrix calculation used in 
Chapter III is repeated here as: (repeated from Eq. 3.86) 

X - M = (« + Tc)(X - p)($ + rc) T 


is not more applicable because it assumes an optimal observer. The co- 
variance matrices M and p of the observer error before and after 
measurement do not represent the modified system, M and p now include the 
additional artificial noise which was applied to the system to obtain a 
better damping of the observer error of the bending mode poles. Thus 
the rms response will be calculated by the following method: 

The state equations of the plant will be augmented by the equa- 
tions of the compensator system to a 2n x 2n open loop repre- 
sentation, The external noise, combined with the measurement noise, 
act simultaneously on the enlarged system. The equations of the 
augmented system are, first, the controlled plant: 


4 x + r Cx. + r w 

Pi pi 2 i 


N(0, Q d ) 


and the observer f 

x i+1 = KHG^ +' [(i-KH) (<HTC) + KHTpC ]x^ + Kv^ v ± N(o,r) 


These equations were derived in Section VI-B, 


(6.18) 


The steady state covariance matrix X of the augmented states is 

a 

the solution of 


X 

a 


T 

$ X $ 
a a a 



T 

a ’ 


(6.19) 
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where 


<T 

a 


* 

P 


r c 

p 


j. 


KM 


P 


(l-KH)(*+rC)+KHT C 


( 6 . 20 ) 


Q 

a 



r 

a 



( 6 . 21 ) 


( 6 . 22 ) 


The observer error gains K, and the state variable feedback gains 

C were calculated in Section VI-C. $ and T include the uncer- 

P P 

tainty for the short period example, and the variation on the bend- 
ing mode frequency. The rms response of the plat states are the 

square roots of the first n terms on the main diagonal of X . 

a 


The results of the foregoing method for the rms response of the 

~6 

modified system (Ax^ = 10 , = 200) is compared to the original 

optimal design, Fig. VI-13. There is a certain increase of the rms 
response of the states. 


RELATIONSHIP BETWEEN THE COMPENSATOR POLES AND SENSITIVITY 

In this section we will generalize the results of Section VI-A. 

It will be shown that there is a definite relationship between the location 
of the compensator^ poles and the sensitivity of the closed loop system. 

The compensator poles are not necessarily well damped or even stable 
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FIG. VI- 13 A COMPARISON OF THE RMS RESPONSE BETWEEN THE 

OPTIMAL DESIGN AND THE DESENSITIZED DESIGN. 

FLIGHT CONDITION No. 8 . Wind gust: cf = 12 

ft/sec: T = 0.5 sec. 

' w 
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for a stable closed loop system. However, our analysis will show that 
if the coefficients of the plant are perturbed, the closed loop pies 
have the tendency to migrate toward the compensator poles. Therefore, 
if the compensator poles are stable and well damped, the closed loop 
system is less sensitive to variations of the plant's parameters. This 
interesting property of the perturbed closed loop system will be proven 
as follows: The characteristic equation of the closed loop system 

(from Section VI-2, Eq, 6,15) is 


«-Iz r c 

p 

KM ( I-KH)($+ t 'c)h-KHT C-Iz 
P P 


o . 


By using (6.13) we will reformulate (6,15), 


x. , = 0 x. + r Cx. 


i+l 


pi pi 


(6.23) 


i+l 


= Ax. + Ky. 

l i+l 


where 


y = Hx 
i+l i+l 


and 


a = (i - kh)(® + rc) 


which yields a different formulation of the characteristic equation 


$ -iz r c 
p p 


KHz A-Iz 


= 0 


(6.24) 


Using the well known relationship, 


det 


A A 
1 2 


A A . 
3 4 


= det A det [A - A A 1 A ] 
4 L 1 2 4 3 J 
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Equation (6,24) will obtain the following form: 


P (z) det[(0 - Iz) - r C 5(^1 KHz] = 0 

c w P P P c (z) 


(6.25) 


where P (z) = |a - Iz I are the poles of the compensator. B(z) are 
C -1 

zeros of (A - Iz) . Let $ = 0 + AO, T = T. The characteristic equa- 

P P 

tion will be 


ju 


B( z) 


P c (z) det j ( 0 + AO - Iz) - TC - - , ^ 

C 


KHz 


] ■ 0 


(6.26) 


or 


P (z) det (0 - Iz) - TC 

c _ PJz) 


KHz + AO 




(6.27) 


Assuming 


AO 


AO . . 

ij 


Equation (6.27) will become 

P c (z)|det[(0-I z ) - r 


rff} det 


[M u ] } 


U\zJ 


(6,28) 


where M. . is the cofactor of M. Now the characteristic equation obtains 
the following form: 


P c (z){det M(z) + det M^} - 0 , 


(6,29) 


However, P^(z)det M(z) in (6,29) is the characteristic equation of the 
nominal system: 
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(6.30) 


P (z) det M(z) = | (<f+rc-lz) [4>(l-KH)-Iz]| £ NP(z). 


All zeros of the NP(z) (nominal poles) are inside the unit circle. 

Finally, the characteristic equation of the perturbed system is 

NP( z) + ..P (z)|m(z)| = 0. (6.31) 

' ij c ' ■ ij' 1 

Recall that P^(z) is the characteristic equation of the compensator. 
Equation (6.31) is now in a familiar form for a root locus analysis, 
where the perturbed coefficient A replaces the gain. An increase 
in A£. will move part of the closed loop poles [zeros of NP(z)] 
toward the poles of the compensator [zeros of P^z)]. Therefore, even 
if the observer error system is stable, it is important to check the 
location of the compensator poles. Unstable compensator poles may cause 
a greater sensitivity to variations in the parameters. 


£. HIDDEN INSTABILITY DUE TO SAMPLING 

In the previous sections we have defined the bending mode as an 
undesirable frequency. It has also been shown that the optimal compensator 
generates a notch filter which attenuates the unwanted frequency. Class- 
ical approaches, such as those given in Borow [BO-1 ] and Sutton [SU-1], 
design the notch filter directly into the feedback loop. The notch 
filter implementation requires mechanizing a second order system on the 
digital computer. Therefore, it seems promising to circumvent the 
necessity of using this filter by taking advantage of the filtering 
properties of the sampler. McGough in 1973 proposed sampling at 

the same rate as the bending frequency and he proved doing this would 
filter the unwanted frequency. The purpose of this section is to show 
that, for certain cases, the closed loop will diverge. This effect cannot 
be detected by the regular stability analysis used by McGough |McG-l ]• 
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In Section VI— F-l we will show what the conditions are to generate the 
instability effect, and in VI-F-2 we will demonstrate it for the short 
period example. 


F-l Problem Statement and Theoretical Background 


The transfer function of the ZOH is given as [RA-lj 


ZOH 


1 


-sT 

- e 

sT 


The frequency response of the zero order hold is described in Fig. VI-14. 


1-e 


- jcoT 


jcoT 



FIG. VI- 14 FREQUENCY RESPONSE OF THE ZERO ORDER 

HOLD (ZOH) 


Intuitively it is obvious that a signal with a frequency near the sampling 
rate or near the multiples of the sampling rates will pass the ZOH 
as a slow varying bias, or d-c. See Fig. VI-15. This idea is especially 
attractive for our principal example where the second bending mode is an 
integral multiple of the first bending mode (25 H, 50 H). Therefore, 


- 116 - 



AMPLITUDE 



FIG. VI— 15 FILTERING PROPERTIES OF THE ZERO ORDER HOLD 


this approach is worth a more detailed analysis. For analyzing this sit' 
uation, we shall use the formulation of Section VI-A-1. The system, 


x 

a 






0 



(6.32) 


y 


i 


Hx . 


l 


where = _ith bending mode state. 

Using an optimal control and no weighting on x^ ... x b^» 
control is 


u . = [C 0 ] 

i a J 


x 


= C x . 
a a . 

l 


(6.33) 
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The exact states x are not known. Any method of reconstruction of the 

a 

states, using $ only, will generate states "contaminated 1 * by the un- 
a 

wanted frequencies 


u. = C x + kx 
i a a. b. ; 

1 i 


k = k(H) • 


(6.34) 


By inspection we can show that the closed loop poles of the subsystems 
are separated. Therefore it will be sufficient to analyze one of 
the second order subsystems. Explicitly written, the subsystem is 


! b - -^bVb - Vb + 


( 6 . 35 ) 


where u is a constant through the interval, T. 

The system (6,35) will be rewritten in a state space form using an 

artificial state x „ This is done in order to close the loop with an 

b 

accelerometer feedback 


V 


r o 




X b 






b b 
2 


0 

O 

- 2 ?b“b- 




u + 


0 


o 


LgJ 


U . 


(6.36) 


Using a ZOH for u we get 


X 

i 



= 

A 

ii 


*1 »2 ° 


*3 *4 ° 

*5 *6 *7. 


‘bJ 




u + 
1 


[ u i “ u ± _ 1 ] (6,37) 


where is a linear function of and of the unwanted frequencies. 

If an acceleration or pickoff is used, then: 
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(6.38) 


u 


i 


C x + k(x. + x ) 
a a b a 


where k depends on the particular structure of the compensator. Thus 
the characteristic equation of (6.37) has the following form 


(t> i z 


<t) -z 

4 


y_k 

V 

<|) 7 -z+gk(l-z 1 ) 


= O 


(6.39) 


The bending modes are nearly undamped (£ = 0.01).- To simplify, we will 

assume that t =0. The various terms and y are: 
s b 


cos w T 

D 


0 = cj sin (j T 

6 b b 




sin U) T 
b 

co 


cos co T 
b 


0 “ to. s in co T 

3 b b 


r, = 


_ JL 


(cos (0 T - 1) 

co b 

b 


0 = cos co T 

4 b 


sin co, T 

co b 

b 


*5 = 0 


If the sampling is done in the same frequency as co^, (6,39) reduces to 


(l-z)(l-z) [l~z + gk(l-z 1 )3 = 0 . 


(6,40) 


The roots of (5,40) are z \ 2 = ^ f Z 3 4 ^ 1^1 > * ^ 

subsystem will diverge. The physical interpretation is as follows. 

The modes x and x are constant at the sampling points, i*e. # x - x 

= x^ But if |gk| > 1 , the acceleration will diverge and the closed 

loop system is unstable. This property can be visualized in Figs, VI-16 


\ 
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VI— 17. 



FIG. VI- 16 HIDDEN INSTABILITY OF AN OSCILLATING MODE 

(kg > 1). 


Diagrams similar to Figs. VI-16 and VI-17 can be traced for -1 < kg < 0 

and kg < -1. Note: (a) For all four cases, x and x are constants at 

b b 

the sampling points. (b) Without using the third artificial state x or 

* b 

an equivalent formulation, this behavior could not be detected. 


Hidden Instability of the F-H Short Period Mode Example Due to 


In order to investigate the effect of sampling at the same frequency 
as the bending frequency, we will use the simulation scheme described in 
Appendix B. The model we will use is the extended short period mode 
described in Chapter IV. That model represents the plant. Investigating 
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FIG. VI-17 ASYMPTOTIC STABILITY OF THE OSCILLATING MODE. 

(1kg | < l) 

the claim of McGough [McG-l] regarding the bending mode which generates the 
notch filter, will not be included in the compensator. The closed loop 
system we will simulate is described in Fig, VI-18. 



FIG. VI-18 SIMULATION SCHEME FOR THE HIDDEN INSTABILITY 
INVESTIGATION. 
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where 


o]; 
oj; 

Using this scheme, the influence of the bending mode is eliminated from 
the compensator. However, the measurements y are combined from the 
output of the whole plant. 

The simulation results plotted in Fig. VI-19 and Fig. VI-20 con- 
firm our analysis of Section VI-F-1. The behavior of the bending mode 
for nominal accelerometer pickup is plotted in Fig. VI-19. We can see 
that the bending mode is marginally stable. 

Figure VI-20 plots the behavior of the bending mode for the case 
in which the accelerometer pickup was increased by a factor of 10. 

The system is unstable. The behavior of the bending mode and the input 
in Fig. VI-20 does not exactly match the schematic behavior charted in 
Fig. VI-17. The reason for this difference is that in Fig. VI-17 we 
plotted the schematic behavior of the homogeneous solution of the 
bending mode. In Fig. VI-20, the rigid body mode is also excited. 



G. SUMMARY 

1. An optimal compensator of a closed loop system having an un- 
wanted frequency is essentially a notch filter. The depth of the notch 
depends on the amount of pickup of the unwanted frequency. 
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2* In the F-H short period example, the closed loop control system 
is sensitive to a bending frequency variation. The sensitivity increases 
for lower sampling rates. 

3 # The sensitivity to a bending frequency variation can be reduced 
by increasing the damping of the observer error poles, 

4, Systems with more stable compensators are less sensitive to 
parameter variations. 

5, Sampling at the bending mode frequency may cause an instability 
which contradicts the results given in McGough [McG-l]. The reason for 
the instability is that McGough didn* t consider the intersample behavior, 

6, The important parameters which cause the sensitivity are poorly 
damped, almost undisturbable modes, and low sampling rates* 
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VII, SAMPLING TIME AND ROUGHNESS OF CONTROL 


Most of the digital systems which control analog (continuous) plants 
use the zero order hold (ZOH) as a reconstruction hold. The abrupt action 
of the ZOH for higher sampling rates is lessened and is further smoothed 
by the filtering properties of the various electro-mechanical actuators. 

The tendency is to shorten the actuator time constants as much as possible 
in order to satisfy various time response criteria. Therefore the action 
of the controls becomes more abrupt at lower sampling rates. However, 
there is no documented analytical approach to control roughness in the 
digital control literature. The purpose of this chapter is to formulate 
and analyze a criterion which will enable us to compare the roughness of 
control for different sampling rates and for different control laws. 

The first experimental evidence of the necessity for such a rough- 
ness criteria appeared recently [MA-1], In that paper, J. Mathews 
describes the Saab digital flight control system. This is a highly 
successful digital control impelmentation. The first flight test was 
made in March 1973, Upon installation of the system in the aircraft, it 
was discovered that the servo valves responded to the small steps in the 
output wave form. It was originally believed that the high sampling rates 
(40 cps in roll and 80 cps in pitch) would be above the bandpass of the 
servo valve. This turned out not to be the case and a lag had to be 
added to the servo preamp. The roll axis computation rate was changed to 
80 cps to minimize the size of the lag filter. This example clearly 
shows the need for an analytical tool which will help the designer estimate 
roughness of control. 

The basic concept in this chapter are the definitions of roughness 
functions (rf) ♦ The fundamental RF will be defined as the weighted sum of 
the squares of the abrupt changes in the states r derivatives or the con- 
trol inputs. 

In Section VII-A we will define the RF for an impulse response of 
a continuous system controlled by a digital controller, A method based 
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on eigenvector decomposition will be given for computing the RF, A 
simple first order example will illustrate the new concept. The RF of the 
F-H short period mode for different sampling rates will be investigated in 
detail. The surprising result is that if the dynamic properties of the 
closed loop are kept constant, the RF is not a monotonic function of the 
sampling interval . Actually, it will be shown that for this example, 
the RF has a maximum. The minimum RF is obviously zero (for T - 0). 

In Section VII-B we will define the mean roughness function, 

(RF^), for a continuous system, disturbed by an external noise and 

controlled by a digital controller. A direct solution of the RF^ wil 

be given. We will also investigate the RF^ of the short period example 

for different dynamic properties of the closed loop while keeping the 

sampling rate constant. An intuitively predicted phenomenon, that the 

RF increases for larger mean values of the control, will be verified. In 
m 

Section VII-C we propose and investigate different RFs which enable us to 
compare roughness of systems controlled via different reconstruction holds 
and actuators. 


A. THE ROUGHNESS FUNCTION /fa 1 , 


A-l Definition 

A continuous linear system, controlled discretely, has the follow- 
ing form: 

x = Fx + Gu 

where u = constant, t 4 £ t < t. + T. Written in a discrete form, 

* x 1 

X = $x + Tu . (7.2) 

l+l i i 

When a step function is applied to the system (7.1) at time t^ f some 
of the variables x abruptly change their magnitude from to x + ^ . 

These are all the variables (not states), which are directly influenced by 
the input u^. In other words, the abruptly changed variables are the 
derivatives of the states for which the coritrol distribution matrix G 
has nonzero elements. We shall limit our discussion to the case in 
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which the state variables x are not abruptly changed; i.e., x + ^ = 
x_^ = x^. This is the common case in mechanical systems, where a step 
input in forces changes the acceleration. The requirement that none of 
the states be abruptly changed is necessary in order to obtain an explicit 
expression for the RF. 

If a full state variable feedback is implemented, the equations 
(7.1 ) and (7.2) can be combined as the following form: 

i - x_ = GC [(0 + TO) - I]x. . (7.3) 

+ i+l i+1 1 


Proof: The input through the time interval T, 

is 


"i - c *i 


(7.4) 


and 


i+1 


Fx.,. + Gu 
l+l 


i+1 


(7.5) 


i+1 


Fx , + Gu 
1+1 + i + i 


Combining (7.5) and (7.4) 


X 

+ 


i+1 


x 

~i+l 


+ “Iv. 



(7.6) 


By using (7.2) we obtain the desired form (7.3). This concludes the 
proof . 

The quantity |x^ i - | indicates the roughness of control in 

the instant t^. We will define a RF of a closed loop system as a 
scalar non-negative number RF: 


RF 


a 



(7.7) 
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where 



t = N X T 

W S 0 . 


In the definition of (7.7), RF = RF(N). By using (7.3), we may express 
RF as a function of x^: 

N 

rf - V x [(*+rc) - i] T c T GW® + rc) - i]x . < 7 * 8) 

1 ^ 

1=0 

If x is given, the RF is a function of x and N: 

RF = V x (G+rc) 1 [(<r+rc)-i] T cVW[4+rc)-i](*+rc) x q . 

i — o ° (7.9) 


Obviously, the expression (7*9) is rather difficult to calculate* However, 
a simple form for RF can be obtained for the limiting case N There 

are two methods for the calculation of RF^. 

Method 1 : Using z-transforms and Parseval's Theorem: 


r S GC(* i+1 - x.) 

(7.10) 

(z) = GC( Iz - l)(lz - $-rc)x o 

(7.11) 

RF = ^ r T (a)Wr(z 1 )dz . 

(7.12) 


7 

where *jr is the unit circle* 

Method 2 : By using a Liapunov function. From (7.2) and (7.8) we obtain 

the following set of equations: 
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(o + rc)x i 


(7.13) 


where 


i+1 

V" T t 

RF = > x.Rx. = Tr(x x P ) 

Z- 1 1 o o o 

i=o 


x - given 
o 


r = [(®+rc) - i] T c T G T wGC[(*+rc) - i] 


(7.14) 


(R S 0) 


and P is obtained by a backward sweep from (7.15): 


p . = (® + rc) P (® 4- re) *t R 

j j+i 


Pj^ — R; j — N - 1 y • ♦ • y 0 # 


(7.15) 


For N the solution for the RF reduces to a simple form: 


RF 


2 X i Rx i = Tr (Vo p) > N_>co 
i=o 


(7.16) 


where P is the solution of 


P 


(« + rc) T p(<r + rc) + r . 


(7.17) 


A- 2 Numerical Solution of the Roughness Function 

The easiest way to obtain the limiting case RF^^^ is to use 
Method 2, which involves the solution of the linear matrix equation 
(7,17). This type of equation is extensively treated in Chapter II 
(Section B— 4). The solution we proposed there was to use the eigen- 
vector decomposition algorithm. This algorithm, directly applicable 


— 129 - 



for (7.17) solves for P. After that, the RF is immediately obtained 
from (7.16). 

A-3 Example of Roughness Function 

We will demonstrate the RF on a simple first order continuous system 
controlled discretely. The purpose of this example is to illustrate 
the new concept. A more realistic example will he worked out in 
Sections VII— A— 3 and VII— B-l. 

The first order system is 


x = -ax + gu 


i+1 

r 


-aT 
e x . 

+ ru. 

i 


(1 - e‘ 

•aT\ g 

a 

cx . . 



(7.18) 


Equation (7,18) combined together yields 


X i+1 


- [•"* ♦ (1 - ."> fK 


a. i9) 


The behavior of the controlled system of (7.19) is schematically des- 
cribed in Fig. VI I— 1 . 

A-4 Roughness Function of the Controlled Short Period Mode for 
Different Sampling Rates 

We will investigate the behavior of the RF for the short period 
mode example described in Chapter III, The open loop poles of the 
short period mode are located at s = -2.6 +_ jl3.5 (flight condition 
No. 8). The closed loop poles have been relocated to s = -15 £ jlO. 
The simplified equation of motion is 
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FIG. VII-1 SCHEMATIC DESCRIPTION OF THE ROUGHNESS FUNCTION 


• 


- - 

q 


q 

« 

a 

= F 

a 

- - 




where u is a discrete input, u^(t) = Cx^, t < T < t^ + ^; mini 

mizes the cost function, J. 



(7.21) 


In our RF calculation we are assuming a perfect estimation of the states. 
This is assumed in order to make a parametric study. In this section 
we will investigate how the roughness changes at different sampling 
rates, while the dynamic properties are kept nearly constant. In 
the next section, we will keep the sampling rate constant but we will 
vary the dynamic properties. 

The dynamic properties are characterized by the location of the 
closed loop poles on the z-plane. In Table VII-1 the closed loop 
poles of the short period example are summarized. In order to provide 
a better engineering insight, the equivalent s-plane poles are also 
listed. Furthermore, the trace of the poles on the z-plane is given 
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Table VII-1 


POLE LOCATIONS OF THE SHORT PERIOD EXAMPLE FOR DIFFERENT 

SAMPLING RATES 



Pole Location 

Control Gains 

Sampling 
Interval 
(sec ) 

z— plane 

Equivalent 

s-plane 

c c 

q a 

Open Loop 

— 

-2.6 ± jl3.5 

0 0 

T = 0 

— 

-14.5 ± J59.6 

0.3 0.6 

T = 0.02 

0.73 ± J0.14 

-15.0 ± jll 

0,273 -0.0101 

T = 0.03 

0.63 ± j0, 68 

-15.5 ± jll. 5 

0.245 -0.25 

T = 0.05 

0.11 ± j0.22 

-16.0 ± jl2 

0.198 -0,603 

T == 0.08 

0.18 ± j0.2 

-16.4 ± j 12.4 

0.143 -0.93 

i-3 

II 

O 

• 

>— * 

0.07 ± jO. 15 

-17.0 ± jl2.5 

0.11 -1.06 

T = 1.2 

0,003 

-0.06 

-17.5 ± J12.7 

0.086 -1.1 


in Fig. VII-2. All the z-plane poles and the feedback gains were cal- 
culated by the discrete synthesis program, keeping the weights in 
the quadratic cost function constant. The continuous case, T = 0, 
was calculated by a continuous synthesis program by Bryson et al [Bft-2 ], 
As seen in Table VII— 1 , the closed loop poles for different sampling 
interval were not significantly changed. The reason for the drastic 
change in the control gains will be explained later. The RF for this 
example was calculated for two different initial conditions: (a) 

q o f °» a o = 0; (b) q = 0, « o £ 0. 

For simplicity we assumed that the weigting matrix W has the 
simple form 


1 


W = 


0 



(7.22) 
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Im 



FIG. VI 1-2 THE CLOSED LOOP POLES OF THE SHORT 

PERIOD EXAMPLE 


This is not a serious restriction as it can be easily proven that for a 
scalar input system, the time history of the roughness (Ax\ ) is similar 
to all the states (influenced directly by the input). Therefore, we 
calculated the abrupt changes in q only. The RF of the short period 
example vs different sampling intervals is plotted in Fig. VII-3. In 
order to make the various RF plots independent of the initial values, 
the RF was ^normalized by the initial conditions as follows: for 




(7.23) 


and for q = 0 
o 
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(7.24) 


( rf ) 2 = • 
a i 

o 

Note: (a) For different weighting matrices, W, the plots of RF in 

Fig, VI 1-3 will have exactly the same form and only the scale of the 
vertical axis will be multiplied by some constant, (b) Although we are 
dealing with a linear system, the superposition principle does not hold 
for the RF, i.e., the RF for combined initial conditions is a quadratic 
form of the initial conditions. 

Without a proper simulation, it will be difficult to explain the 
behavior of the RF for different sampling intervals. Therefore a 
complete digital simulation has been done and the simulation results 
are plotted in Figs VII-4 through VII-7. The simulation scheme is 
described in Appendix B, The different plots are organized by the 
following order: 

T = 0,02 sec 

\ Fig. VI 1-4 
T = 0,05 sec J 

T = 0.08 sec > 

> Fig. VII— 5 
T =0,12 sec J 




By inspecting these plots we can see the following properties. 

(a) For 0^=0 (Figs. VII-4, 5), as the sampling interval T 
increases, the initial input u^ decreases. This behavior 
stems from the fact that the optimal gain c decreases for 
larger sampling intervals (see Table VI-1). That change in 
c is easily explained by using the fundamentals of the optimal 
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control theory: the optimal control gain minimizes the quad- 

ratic cost function. By keeping the weights in the quadratic 
cost function constant, the squared area under the u /q q curve 
should be approximately constant. For longer sampling inter- 
vals, the first input acts on the system longer; therefore, 
its amplitude has to be smaller. 

(b) The optimal design essentially generates a damping augmenta- 
tion (see Table VII-1). Thus if q - 0 and (X = 0, the input 

o 

builds slowly until q reaches maximum (Fig. VI-6). But 
for longer sampling intervals, the optimal control doesn’t 
have the time to wait for larger q. Therefore it predicts 
the larger q and generates feedback by increasing (see 

Table VII-1 and Fig. VII-7). 

(c) The fact that the RF is lower for sampling intervals and 



FIG. VII- 3 THE ROUGHNESS FUNCTION OF THE SHORT PERIOD 

EXAMPLE VS DIFFERENT SAMPLING INTERVALS. 
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FIG, VII-7 


TIME 


0,3X10 (sec ) 


T = 0,08 sec 



0,4 


t (sec ) 


0,2 


0,3 



) = 0.5X10 4 (sec” 4 ) 


T = 0.12 sec 


longer than 0,1 sec does not mean that we are going to recommend 
the use of this low sampling rate (Fig, VII-3), There are 
other considerations in the selection of sampling rate, e,g a , 
the time response. 


B, THE MEAN ROUGHNESS FUNCTION OF A CLOSED LOOP SYSTEM DIS- 

TURBED BY AN EXTERNAL NOISE ~~ 


We are also interested in estimating the average roughness of a 
continuous system controlled discretely and disturbed by a random noise 


B— 1 Definition of RF 

m 

The disturbed system is 


= Fx + Gu + w 


w N(0, Q) 


x 


i+1 


<T x . + w. 

Cl 1 


w i -+ N(o, Q d ) . 


(7.25) 


* is the closed loop transition matrix. Q is the covariance matrix of 
c d 

the discretized white noise with a power spectral density matrix Q. 

The mean roughness function, RF^ , of this system will be defined 
as 

RF = Tr E (x, - x . ) (S - x ) T W ) i -* « 

“ * i _1 + i - 1 ' (7.26) 


where Wa 0. Using relation (7.6), we have (repeat) 

Ai = x - K - G ( u i+ i ' V * 

i+l i+1 


Therefore , 

' *i )5 • 

Combining (7.25) and (7.27, we get 


(7.27) 
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E {Ax i+ i} = E{GC($ c X. - x. + Wi )}, w. ^N(0, Q d ) 

where E is a linear operator. Furthermore, x^ and are uncorrelated 

by definition. Therefore 


e{ax. ,} = e{gc(<& - i) x . j 

1+1 C 1 


(7.28) 


where x i is a GaussrMarkov random process with a covariance matrix 

X - Ax, is a linear combination of x. , thus. Ax. is also a Gaussian- 
i i l i 

Markov process with a covariance matrix, X r ^ given by 


X = GC(fl -l)x (® -I)Vg T + GCQ C T G T 
r . c i c d 


(7.29) 


T 

where X. = Efx.x.}, The steady state value of the covariance matrix 
l c i 


X. is the solution of 

l 


X = ® X<3> + <3 . 

c c d 


(7.30) 


(see Chapter II). By combining (7.26) and (7.29) we have the final 


expression for RF 


m 


RF =s Tr(x W) 
m % r ' 


(7.31) 


where 


X r = GC(^ c -I)x(<I c -I) T C T G T + GCQ d C T G T . (7.32) 


Similar to the RF for an impulse response, the RF^ involves the solution 
of a matrix linear equation, (7.30). 


B-2 RF of the Controlled Short Period Mode 
— m 

In this section we will investigate the RF for the short period 

m 

example described in Chapter IV. The only difference between this model 
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and the one used in Section VI-A-3 is in the augmentation of the states 
by the wind model as we will now describe. 

The system: 


+ T u + r 2 w, w -> N(0, Q) . (7.33) 


The open loop poles were located at s = —2.6+_ 1 jl3.5 and the nominal 

weights (Ch. VT-C) relocated the closed loop poles to s = 10 _+ 15. As 

in Section VII-A-3, we are also assuming perfect measurements of the states. 

This is done for simplicity of calculation as we are interested in the 

behavior of the RF for different weights (A A , ) rather than in 

m a/ £> Oif 

its absolute value. The wind gusts model used in this example is the model 

of flight condition No. 8 (o = 12 ft/sec t = 0.5). All the calcula- 

N w 

tions were made for one sampling rate (T = 0.05 sec). 

By varying the values of A ^ and A^y^ we will consider the 
behavior of the following quantities: (1) q^, covariance of the state q; 

—2 

(2) ctj,, covariance of the state 0^,; (3) u f covariance of the control 

u. Furthermore, the weighting matrix W in the RF^ criteria (6.26) 
will have the form 




- 

- 


- 

q 



r 

w 


q 


= 




“T 

* 

W 

L gJ 


0 

« 

w _ 


w 

L e J 


W 


10 0 
0 10 . 
0 0 0 . 


(7.34) 


and we will consider the two components of RF : (repeat of Eq. 7.31) 

m 


RF = T (X W) 
m r r 


namely, RF f 

q 


component of RF , corresponding to 

m 


W 2 ; 
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2 

component of RF corresponding to (Z$) , The RF calculations are 
m m 

as follows: 

The closed loop pole location as a function of different 

weight is plotted in Fig. VII-8. The nominal design was chosen 

in Ch. IV. The behavior of the above listed quantities as a 

function of A / is plotted in Fig. VII-9. The covariance 
q/B 

of q (q2) decreases for larger weights A # • But this cannot 

q/ B 

be done without a penalty, which (a) an increase in RF , 

_2 ^ 

(b) an increase in u . 


The behavior of the system as a function of A^^ is more 
complicated. As seen in Fig. VII-10, only q^ increases as 
increases. The physical explanation is given in Fig. 
VII-8 where the increase in A^ 0 keeps nearly the same damp- 
ing but decreases the natural frequency; this means that the 
systems reponds slowly and there are less abrupt actions in the 
control. As a conclusion of this section, we may say that the 
RF^ provides us with additional insight into the system. 


A a/B 


C. ROUGHNESS FUNCTIONS FOR FIRST ORDER HOLD, TRIANGULAR 
HOLD, AND PLANTS WITH FILTERED INPUT 

The object of this section will be to suggest and investigate dif- 
ferent functional forms which will enable us to compare roughness of 
control for different reconstruction holds and for plants which include 
actuators or filtered inputs. 

The properties we expect from a roughness function are: 

(a) The RF should be related to the difference of numbers in 
the discrete input sequence; 

(b) For T = 0, the RF is zero; 

(c) In order to establish experimentally a relationship between 
the RF and the pilot T s evaluation, the RF should be repre- 
sented by a positive function; 
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FIG. VI 1-8 SHORT PERIOD MODE. Pole locations as a function of weights 
on q and Flight Condition No. 8. T = 0,05 sec. 







(d) The RF has to be user-oriented and it should have a conven- 
ient computer program form. 

We will investigate the RF for four different control configura- 
tions: (a) zero order hold/zOH; ( b) first order hold/FOH; (c) tri- 

angular hold; (d) filtered input (or actuator lag)* The block diagrams 
of the system and their schematical time responses are plotted in Figs* 
VII-11 and VII-12* It should be noted that in the classical sample data 



FIG* VI I- 11 DIFFERENT RECONSTRUCTION HOLD CONFIGURATIONS 
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ZERO ORDER HOLD 



T 2T t 



FIG. VII-12 THE SCHEMATIC BEHAVIOR OF DIFFERENT RECONSTRUCTION HOLDS 
AND THE PARAMETER DEFINITIONS FOR THE RF CALCULATIONS. 
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theory, the triangular hold is an unrealizable reconstruction hold and is 
used only for non-real time integration. However, if optimal control 
is used, the input u is actually a part of the state vector and is 
fed back (see Fig. VII-ll). Therefore the triangular hold is realizable 
in this formulation. 

The proposed RFs are: 



These four RFs can be calculated by simulation only. RF and RF are 

1 2 

actually identical to the RF which was defined in Section A for a 

zero order hold. Only for this case does RF and RF have an analytical 

i. z 

solution , 

5* COMPARISON OF ROUGHNESS OF THE F-H MODEL FOR DIFFERENT 

RECONSTRUCTION HOLDS 

Each of the four RFs will be demonstrated for the F-H short period 
mode (Flight Condition No. 8, T = 0.05s) using different reconstruction 
holds and actuator. The optimal feedback gains and the digital simula- 
tion was done by the Program DISC and by the simulation program described 
in Appendices A and B. Since the two programs are constructed for a 
ZOH (zero order hold), the FOH (first order hold)— the triangular hold-- 
and the filtered input were reformulated to a ZOH equation: (repeat) 

x = Fx + Gu 
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where u is constant for 0 ^ t < T. For the FQH 


r~ ~i 


- 

— 




- 



■* *" 

X 


F 

G 


X 

+ 

G 

0 


* 

U 2 

u f 


0 

0 


u 


0 

I 


k 

j x 


— 

- 


- lu 


- 



- - 


where 

d 


u = 


constant 1 
constant 
6 , + 


0 ^ T < T 


The triangular hold; 


- - 


- 

- 


- - 


- - 

X 


F 

G 


X 


0 


= 







u 


0 

0 


u 


I 

— - 


_ 

- 


^ - 


^ - 


where 6 = constant, 0 ^ T ^ T. 
The filtered input; 



(7.35) 


(7.36) 


where u = constant, 0 ^ t ^ T. 

The feedback gains (weighting matrices) were so chosen that the closed 
loop poles of the short period mode were placed in identical positions 
for all the four control configurations (z = 0.44 + j0.22). 

To obtain a visual impression of the roughness, the behavior of the 
input u is plotted in Fig. VII-13 (for an initial condition in a). 
However, from visual inspection, we cannot definitely conclude which hold 
is ’better 1 . The pilot’s rating is necessary to evaluate the roughness 
and to relate it to the proposed roughness functions. The numerical 
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0.1 0.2 0.3 


FIG. VII-13 TIME HISTORY OF THE CONTROL INPUT OF THE DIFFERENT 
RECONSTRUCTION HOLDS OF THE F-H, Flight Condition 
No. 8 , T = 0 .05s 
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values of the proposed roughness functions are summarized in Table VI 1-2. 


Table VII-2 

THE ROUGHNESS FUNCTIONS VALUES FOR THE F-H SHORT PERIOD 
MODE OF THE DIFFERENT RECONSTRUCTION HOLDS 



ZOH 

FOH 

Triangular 

Hold 

Filtered 

Input 

H 

1.42 

0.95 

0.16 

1.1 

rf 2 

1.42 

0.9 

0.36 

1.0 

CO 

Pn 

0.9 

0.6 

0.3 

0.9 

RF 4 

0.9 

0.62 

0.6 

0.8 


From Table VII-2 we may conclude that for the four proposed RFs, the 
triangular hold has the lowest value. 


E. SUMMARY 

1* The experimental results of Mathew for the digital autopilot 
of the Saab Wiggen airplane [MA-l] indicates the need for some measure 
of roughness as a design tool. 

2. We have defined such a measure, calling it the "roughness func- 
tion" (RF), which can be calculated by relatively simple algorithms. 


-152- 







VIII. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


A-l Selection of Sample Rate 

Four principal factors influence the selection of the sample rate 
for the discrete control of a continuous system. These factors are: 

a. the time response to control inputs, 

b. the response to an external noise, 

c. the sensitivity to- variations of parameters, 

d. the roughness of the response to control inputs . 

Each of these factors and its relation to the sampling rate should be 
carefully investigated. 

a. The time response . 

The time response is directly related to the sampling in- 
terval. For future fly-by-wire controlled aircrafts, there is an in- 
herent time lag between the pilot’s decision and the execution of his 
command. This was discussed in Ch. III. 

b . Response to an external noise . 

The response to an external random disturbance is a function 
of the sampling rate. It was found that in discrete systems, there is 
a limit on noise alleviation. This limit is a function of the sampling 
rate. 
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c. Sensitivity to variations of parameters 

The sensitivity to variation of parameters can be reduced by 

a modification of the discrete compensator. In the F-H short period example 

the bending mode frequency is one of the major uncertainties. It was found 

that the optimal compensator generates a narrow and sensitive notch filter 

which attenuates the bending frequency. For a 10 percent error in the 

bending frequency, the system is unstable for lower sampling rate (for 

w /w < 2). By modifying the compensator to broaden the notch, the bend- 
s b 

ing mode was stabilized for much larger errors in bending frequency. 

d. The roughness of response to control inputs 

The roughness of response to control inputs is an inherent 
property of the zero order reconstruction hold. The roughness function 
defined in Ch. VI is used to give an analytical measure of the rough 
behavior of the system. This criterion, demonstrated on the F-H short 
period mode, will help designers to evaluate the behavior of a controlled 
system and select a suitable sampling rate* 

A- 2 Relationship Between Sensitivity an d Compensator Poles 

As is already known; the stability of a closed loop system is 
unrelated to the stability of the compensator. A stable system may have 
an unstable compensator but it has been found that if the compensator 
poles are more damped, the closed loop system has lower sensitivity to 
variation of parameters. 


A— 3 Filtering Property of the Zero Order Hold 

The zero order hold has a well-known property: it filters all 

the frequencies which are integral multiples of the sampling rate. 
However, it has been found that if the system includes a subsystem 
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oscillating in the same frequency as the sampling rate, the whole system 
may be unstable in the intersample intervals. This instability cannot be 
detected by a discrete analysis. It is necessary to investigate the 
intersample behavior. 


B. RECOMMENDATIONS 


Further research is required on the following topics. 

1. A theory and a computer algorithm for synthesis of an optimal 
discrete system which includes a sensitivity minimization. 

2. An analysis of the behavior of a closed loop system which 
uses the modified first order reconstruction hold (proposed in Ch. II-A). 

3. An experimental investigation of the relationship between both 
handling qualities and pilot comfort, and the numerical values of the 
roughness function. 

4. An analysis of the reduction of roughness by filtering the re- 
construction hold output, and of the resulting influence on the stability 
and complexity of the system. 

5. A method which will help the designer decide on the minimal 
representation model of a given system. This representation should be 
sufficient to construct a control law for the actual plant. The behavior 
of this controlled plant will be unaffected by the unmodeled states. 
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Appendix A 


PROGRAM DISC : OPTIMAL DISCRETE CONTROL AND FILTER 

SYNTHESIS BY EIGENVECTOR DECOMPOSITION 


In this Program there are various options. The basic calculations 
are as follows. (a) For a system given as 


x._ = + r u , 

i+l i 1 i* 


(A-l) 


DISC calculates gain C. u. = Cx. which minimizes the cost function 

7 i l 


T 1 v* T a T^ 

J = 2 L x i Ax 1 + U.BU. 


i=o 


A > O 
B > 0 


( A-2) 


( b) For a system given as 


x. (1 = Qx. + r o w. 

1+1 i 2 1 


y. = Hx. + v. 

i li 


w i N (°> Q d ) 

V ± -* N(0, R) , 


(A-3) 


DISC calculates the matrices M, P, K, where K represents the 
gains in the Kalman filter. 


x. , , = Qx. 

i+l l 


(A-4) 


= x ± + K( y i - Hx ± ) 


and M = the error covariance matrix before the measurements, and 
P = the error covariance matrix after the measurements. 

(c) For a system combined from controller (a) and filter ( b) , DISC 

SECEDING PAGE BLANK NOT FILMED ~ l57 ~ 




calculates the covariance matrix X of the states, and the covari- 
ance matrix U of the control. X and U are the solution of 


x - m = ($ + rc) (x - p) (& + rc)' 


u = cxc . 


Various Options 

The various options are ( a) For 


x = Fx + G^u + G 2 w w N(0, Q) , 


and if the controls are reconstructed by a zero order hold, DISC 


calculates 


$ = e 


i - i « F \< 

•'a 


Q d = ^ 'r(x)G 2 QG^<l(T) T dT 


and sets r 2 = I (unity matrix). T is the sampling interval. 


( b) If a continuous cost function, 


00 

J 1 =: ~ ^ (x T A t x + u T B t u)dt 


is given, DISC calculates the control C, which minimizes the cost 
function J* by transforming J’ to J by 


* “IP 

_ A A ^ x 

t-M " ” 

i-o L J . 

> A 2 2 J [“ 
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(r)A<r(T)dT 


1 T 

= \ <T (r)Ar(T)dr 

J o 

T 

= ^ [r t (t)at( p ) + B]dr 


A = A , 
21 12 ' 


and transforms the system (A-l) into 


X i+1 ^ " r i A 22 A 21' >X i + r i U i 


'1 = U i + A 22 A 21 X i = C ' x l ’ 


After this transformation, (A-9) obtains the form of (A-2), i.e., 


J ^ X i ( ' A ll ' A 12 A 22 A 21' >X i + U i A 22 U i‘ 


DISC solves for C 1 (for A-ll and A-13) and calculates the desired 
C from 

c = c ’ - All 

(c) If perfect measurements are assumed, DISC calculates the covariance 
of the states, X, from 

X = (« + Tc) X (<T + rc) T + Q d . (a-i 


(d) DISC calculates the open loop eigenvalues and eigenvectors, (The 
closed loop eigenvalues and eigenvectors, and the estimate error 
eigenvalues and eigenvectors, are always calculated and printed,) 
Using DISC 
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DATA 


FORMAT 


CARD 7 IOC/ IQ/ INQ/ ICC/ 


Format (412) 


8 

9 


ns/ nc/ nz/ ng/ t 


F or $ 
G 1 or r 
A 
B 

G„ or T 


2 *2 
Q or Q 


-\ 


• H 

/* (last card)^ 


Format (413, F10 # 5) 


Format (6F12.5) 

row by row 
each row starts 
by new line 


CARD 7 : 

IOC 

IQ 

INQ 

ICG 


OPTIONS 

= l open loop eigenvalue and eigenvectors 
~ 1 response for perfect measurements 
= 1 response for estimated measurements 

= 1 continuous cost function 


If one or more of these options are not desired, set values to zero. 


CARD 8 


NS 

number 

of 

states 

NC 

number 

of 

controls 

NZ 

number 

of 

measurements 

NG 

number 

of 

disturbances 

T 

sampling 

interval 
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Important ; In order to simplify the utilization of this program, 
all options are activated by necessary parameters and data. All 
unnecessary information should be deleted. 

The various options and instruction for data input are summarized 
in Table A-l . 


Table A-l 
INPUT TO DISC 



C 

A 

R 

0 

CONTINUOUS SYSTEM 
PLUS 2 OU 

DISCRETE SYSTEM 

REGULATOR ONLY 

7 

■ 

0 

0 

X 

■ 

■ 

B 

BB 


8 

M 

NC 

"o'" 

0 

B 



■B 

0 

9 

B 

3^ A B 

# a u 

REGULATOR PLUS 
RESPONSE TO EX- 
TERNAL DISTURBANCE 

n 

B 


B 

B 

■ 

X 

mm 

0 


H 

g 


m 

g 

■ 

NS 

a 

a 

a 

0 

i 

F G r A B G 2 Q 

* r i A 8 P 2 % 

FILTER 

n 

X 

0 

0 

X 


. 1 

0 1 0 
i 

° 


6 



EM 


D 

Efl 

a | nz 

NG 

0 

9 

1 F C 2 R H 

— 

® r 2 r h 

REGULATOR PLUS 
FILTER 

7 

. 

X 

X 

X 



* 

1 

X 

0 


8 

mm 


ESI 


T 

KBI 

XC ! 

SZ 

NG 

0 

m 

F CjA B G 2 Q R H 

<5 r 1 A E ^ Q d R H . | 


x 0, 1 ... different options; 9^ order of reading matrices in. 

Note : If more than the minimal amount of necessary parameters and 

data are read in, the program execution will fail. For exam- 
ple, you cannot: 


(a) 

set 

ICC = 1. 

for pure discrete system 

(b) 

set 

IQ = 1. 

and INQ — 1. simultaneously 

(c) 

set 

NC ^ 0 

if you are asking for filter 

The following 

is the 

computer 

program for DISC. 
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o# 


E 


F 1 
o 

g 0 

Is 

a* 




to 

I 


6 » 
7. 

8 ® 
?. 

10c 

11c 

12 . 

1?. 

14. 

15. 

16. 
17. 
IB. 
19. 
2C. 
21 . 
22 . 
22.1 
22.2 

22. 3 

22.4 

22. 5 
22.6 

23. 

24. 

25. 

26. 

27. 

28. 

29. 

30. 


//CISC JOB *5270, 202, 2., 10* ,K ATZ,*SGL FVr L*1 
// EXEC FORThCL, UVEL=G 
//FORT. SYS IK DO * 

IMPLICIT P.EAL*0<4-l-tO-Zl 

C THE FOLLOWING ARE DIMENSIONED SC THAI AN 16 ORDER OR SMALLER SYST 
C WITH UP TO 16 CONTROLS CAN FUN 

DIMENSION A<32,?2),RM(32,32),WR(22),WI(32),X<32,22>,PRO(i2,22), 
1CWRU6 ) • Chi (14) , ICNT( 32 ),C( 32) , INT(32) , RA( 16, 16) ,ACL( 16,16) , 

2GN ( 16,16) , QU6 ,16 ),RVEC( 16 ,16) ,VECFN( 16, 1 6 ) , VEC IN ( 16 , 1 6 ) , 

BY EC I ( 16,16) , VRV ( 16 ), V IV C 16) ,VPRV( 1 6 > , V» I V ( 1 6) , W 1 1 ( 1 6 , 1 6 ) , C < 1 6 ) , 
9LT< 256) ,MT (256) ,6( 16,16) , GPU 6, 16 ) , 

3W21 ( 16 , 16) ,TCB( 32, 16) , ATI 256) ,C I ( 16) ,CT (16,16) ,61(16,16) , 
5Y(16,1 6 ) ,FT( 16,16) ,FT IA< 16,16 ) ,FTIY( 16, 16 ), FT IYA (16,16) , 

7SC( 16, 16) , V F C R ( 16,16) , 

40(16,1 6 ) , FBGC (16,16) , VFC ( 16 , 16 ) , V F ( 32 , 3 2 > , H ( 1 6 , 1 6 ) , 

1 AQ( 20, 20) ,W1Q( 20,20 > ,W3G( 20 ,20 ) ,CC (20,20 , 

1 FTQ( 20,20 ),Q21( 16,16) ,A 02 l (16,16 ),GA( 16,1 6) ,00(16, 16), 

1R (16,1 6 >,RI< 16,16) ,BAK< 16,161 ,GL( 16, 16) , GAT (16, 1 6), 

1GAT1 ( 16,16) , G I < 16,16) ,FBGCK(16, 16 ) ,QOGi 20 ,20 
REAL*4 F M T ( 1 6 I 

CALL ERRSFT (207, 256, -1,1,1, 21 5 ) 

C 

C 

C THE 'MAIN ' PROGRAM READS THE DIMENSION , THE OPTION CARO 
C AND CALLS FOP EXECUTION SUBROUTINE. 

C 

c 

21 READ) 5,7432 ,END=19 ) NS, NC ,NZ , NG , T 

NC1=NC 

IF(NC.GT.G) GO TO 1 
NC=NC+ 1 

1 NG1= NG 

IF (NG.GT.O) GO TC( 2 
N G=NG + 1 

2 N Z 1 =N 7 





I 

(-■ 

0 
w 

1 


21 . 

32. 

33. 

34. 

35. 

36. 

37. 

38. 

39. 
<0. 

41. 

42. 

43. 

44. 

45. 

46. 

47. 

48. 

49. 

50. 

51. 

52. 

53. 
53.1 


IF (NZ.GT.C) GO TO 3 
NZ = NZ + l 

3 READ(5,51> IOL, IQ, INQ,ICC 

51 FORMAT (41 2) 

7432 FORMAT (413. FI 0.5) 

N2=2*NS 
MHS=NS#NS 
NNC=NS +NC 

CALL £XEC(NS,NC , N2 ,A, RM , WR f Vi I , X » PRC , CWft ,CWI , ICNTj»D, I NT, 
1ACL,$C,GN,Q,B,BI ,G,FBGC ,GM,BA,C,CI ,CT,MT,LT ,AT,W21 ,TCB,W11,VF, 
3NC1 ,NG1,NZ1, Y,FT,FTIA,FTIY,FTIYA,H,NZ, 

2 VEC, R VEC, VECRN,VECIN,VECR,VECI,VRV,VIV,VRRV ,VRIV,MHS, I A, 

4 A G» Wl C,W3C,CQ,QDC» FTQ ,N NC , Q21 , AQ2 1 ,NG, T ,GA , QO , R , R I , BAK , 

4FBGCK »GL*G4T»GAT1,GI » I0L»I C * I NO. ICO 
GO TO 21 
10 RETURN 
END 

SUBROUTINE EXEC ( NS ,NC t N2 » A , RM , WR , Vi I , X ,P RO , C WR , CW I , I CNT , D , I NT, 
1ACL»SC»GN»Q»B,BI ,G,FBGC ,GM , BA ,C ,C I ,CT , MT , LT ,AT ,W2 1 ,TCB,W11,VF, 
3NC1 ,NG1 ,NZ1,Y,FT,FTI A,FTIY, FTIYA,H,NZ, 

2 VEC, R VEC, VECRN, VEC IN, VECR,VEC I , VR V ,VI V, VRRV ,VR IV ,MHS, I A, 
4AC,WlC,k3Q,CQ,QDG,FTQ,NNC,Q21, AQ2 1 ,NG, T ,G A , CD , R , R I , BAK , 

4FBGCK ,GL, GAT ,GA T1 , GI , IOL ,1 C ,INO , I CC) 


53.2 
53. 3 
53. 4 

53. 5 

54. 

55. 

56. 

57. 

58. 

59. 

60. 
61. 


C 

C •EXEC" READS ALL THE MATRICES FIRST... 

C 

C 

IMPLICIT REALMS < A-H,G-Z > 

DIMENSION A(N2»N2) ,RM(N2,N2),WR(N2)»WI(N2) ,X(N2,N2) ,PR0(N2»N2 » , 
lCWR(NS) ,CWI (NS) * I CNT ( N2 ) » D( N2 ) » INT (N2 ) » BA ( NS »NS ) , ACL( NS, NS ) * 

2GN( NS ,NS> , Q(NG ,NG ),RVEC(NS ,NS> , VECRN( NS , NS > , VEC IN ( NS, NS » ,VECP(NS, 
7 NS ) , 

8 VEC I (NS, NS) , VRV (NS) , VI V ( NS ) ,VRRV ( NS) , VR I V ( NS) , Wi i (NS ,NS > ♦ C ( NS > , 
9LT( MHS ) ,MT( MHS) , B ( NC , NC ) , GM (NS , NS ) , 

3V»21(NS ,NS > ,TCB( N2 ,NS) , AT(MhS) ,CI(NSI ,CT (NS , NS ) ,61 <NC,NC) . 



I-* 

2 

i 


62. 

63. 

64. 

65. 

66 . 

67. 

68 . 

69. 

70. 

71. 

72. 

73. 

74. 

75. 

76. 

77. 

78. 

79. 
GO. 
81. 
82. 
82. 
84. 

65. 

66 . 

87. 

88 . 
86.1 
88.2 

89. 

90, 

si. 

92. 

S3. 

94 • 

95 • 


5 Y INS, NS), FT! NS, NS) ,FT IA l NS , NS ) , FT I Y ( NS , N$ ) , FT I YA ! NS , NS ) , H I NZ , NS I , 
4G!NS,NC > ,FBGC!NC,NS) , VFC!NS ,NS ) ,VF !N2,N2 I , SCINS,NS I , 

1 AQ( KNC, NNC) ,W1Q ! NNC, NNC > , ^3G( NNC, NNC) ,C 0! NNC ,NNC ) , 

IF TO (NNC, NNC) , Q2 1 ( NC,NS ) , AG21C NC ,NS I ,GA( NS t NG) , 00 ( NS , NS I , 
1R(NZ,NZ),RI(NZ,NZ> , BAKINS ,NS> ,Gl( NS, NC> ,GAT(NG,NS>, 

1GAT1! NGfNS ) ,GI! NS ,NS) , F BGCK (NS , N2 ) ,CDQ< NNC » NNC ) 

RFAL*4 FMT( 20) 

C READ ALL MATRICES 

C 

C 

I CN=0 
NC -NC 1 
NG= NG1 
NZ=NZ1 

DO 61 1*1, NS 

fci REACf 5,7433 MBA! I ,J) ,J=1,NS) 

IF (NC.EQ.O) GO TO 10 
DO 62 1=1, NS 

62 READ! 5, 7433 MG( I , J >, J=1 ,NC) 

CALL MAKE I GL,G» NS »NC) 

00 63 1=1, NS 

63 READ(5»7433) (SCI I » J) ,J=1,NS) 

DO 64 I=1,NC 

64 REAC(5, 7433MB! I, J ),J = 1,NC) 

7423 FORMAT! 6G12. 5) 

3 FORMAT!' ','OPEN L GOP ' DYNAM ICS MATRIX....') 

4 FORMAT! 8( 2X, 1PD13.6M 

WRITE! 6, II 

1 FORMAT Cl',' NEW CASE *) 

WRI TE ( 6,6666 ) NS,NC,T 

6666 FORMAT! '1*,* THE ORDER OF THE SYSTEM =» , 12, », NUMBER CF CCNTRQ 
1 L S= * , l?,', SAMPLING TIME = * ,F10. 5) 

IF { ICC.EQ.O ) GO TO 45 
WRI TE (6,111) 

111 FORMAT! 'O', * YOU APE USING THE C0N7. COST FUNCTION ',//) 

45 WRI TE !6 ,3 ) 



-S9I 


5 6 . 

97. 

98. 

99. 
100 . 
101 . 
102 . 

103. 

104. 

105. 

106. 

107. 

108 . 

109. 

110 . 
111 . 
1 12 . 

113. 

114. 
1 15. 
116. 
1 17. 
118. 

119. 

120 . 
121 . 
122 . 

123. 

124. 

125. 

126. 

127. 

128. 
129. 

1 30. 
131 . 


00 5001 NI=l,NS 

5001 W R I T E ( 6 , 4 ) ( B A ( N I , NJ ) , N J= 1 , NS ) 

ViRITE(6,7004) 

7004 FORMAT ! '0 ' , ' THE CONTROL DISTRIBUTION MATRIX IS...',//) 

DO 7005 NI=1,NS 

7C05 W R I T E ( fc » 4 ) ( GIN I , N J ) ,NJ = 1,NC ) 

WRI TE(6, 3041 ) 

3041 FORMAT ( 'O',' THE STATE COST MATRIX IS....',//) 

DO 7007 NI =1 » NS 

7007 WRI TE( 6,4) ( SC I N I , N J) ,N J = 1 , NS ) 

WRITE! 6,3042 ) 

3C42 FORMAT! 'O', 'THE CONTROL COST MATRIX IS. ...»,//) * 

DC 7009 NI = 1 , NC 

7009 WRITE (6, 4) C 8! N I , N J) ,NJ = 1, NC ) 

10 IF (NG.EQ.O) GO TO 15 
DO 72 1=1, NS 

72 REAC! 5, 7433)1 GA( I , J > , J= 1 ,NG ) 

DO 73 I =1 , NG 

73 REA C! 5 , 7433 ) (Q( I , J > , J=1 ,NG ) 

IF INZ.EQ.O) GO TO 11 

DO 74 1=1 t N Z 

74 READ ! 5 , 7433 ) t R( I,J),J=1,NZ) 

DO 65 1=1, NZ 

65 REAC(5»7433)(H( I,J),J=1,NS) 

11 WRI TE (6,6667) NS,NG,NZ,T 

6667 FORMAT! '1',* THE ORDER OF TFE SYSTEM , 12, ', NUMBER OF PPCCES 
IS DIST. = ', 12,', NUM. CF MEASUR. = • , 1 2 , ' , SAMPL. TIME = ',F10.5) 
IF (NC.CT.O) GO TO 12 
WRITF(6,3> 

DO 8001 N 1 = 1 ,NS 

8C01 WR I TE ( 6 ,4 ) ( BA( N l , NJ > , N J= 1 , NS ) 

12 WR ITE<6» 8004) 

3C04 FORMAT! 'O', ■ THE STATE NOISE DISTRIBUTION MATRIX IS....*,//) 

00 8005 NI =1 , NS 

8CC5 WRITE! 6 , 4 ) ( GA ( NI , NJ ) ,NJ=1 ,NG) 

WRI TE (6,8041 ) 



1 32. 

8041 

FORMAT! 'O' ♦ 1 THE STATE NCISF POW. $• DEN. 

133* 


DC 800 7 Nl =1 « NG 

13^. 

8C07 

WEI Tc ( 6 ,4 J <Q(NI »NJJ,NJ = 1,NG) 

13?. 


IF (NZ.EQ.O) GO TO 15 

1 3fc* 


WR!TE<6,8942 ) 

137. 

3C42 

FORMAT! *0" , 'THE MEASUR. NOISE POW. S. OEN 

128. 


OC 800° NI =1 ,NZ 

139. 

8009 

WRITE! 6, 4) { R(NI, NJ) ,NJ=1»NZ> 

1A0. 


IF (T.EG.OI GO TO 68 

143. 

68 

WRI TE!fc,3043) 

144. 

304? 

FORMAT! 'O' , ' THE MEASURMENT MATRIX F IS... 

145. 


DC 7010 NMItNZ 

146. 

7C10 

WR I TE (6,4) (H(NI»NJ),NJ=1»NS) 

147. 

15 

CONTINUE 

1 48. 


IF <NC.GT.OJ GO TO 6 

149. 


NC-NCl+1 

150. 

6 

IF ING.GT.O) GO TC 7 

151. 


NG=NG1+ 1 

152. 

7 

NCNC=NC*NC 

152. 1 

C 


152. 2 

C 


152.3 

C HERE ! OPTION) EXEC CALCULATES EIGEN VALUES 

152.4 

C 


152. 5 

C 


153. 


IF ( IOL.EQ.O) GO TO 5 01 

1 54, 


DO 53 I = 1 » NS 

1 55. 


DO 53 J = 1,NS 

156. 

53 

GN! I, JJ * BA(I, J) 

157. 

£**************«******** *** *************** 

158. 


CALL BALANC ( NS , NS ,GN , L QW , I F I GH ,0 ) 

159. 


CALL ELMHES INS ,NS .LOW. IF IGF. GN , I NT ) 

160. 


CALL H0R2 ! NS»N$»LCW» IHIGH.GN, CWR.CWI ? A CL 

161. 


CALL ELM8AK < NS , LOW, I H I GH , NS, GN , INT , ACL ) 

162. 


CALL BALBAK (NS ,NS ,LOW , I H I GH, NS , D , ACL ) 

163. 

(^ft*******^*********************** ********* 

164. 


WRITE ( 6 ,2030 ) 


MATRIX IS* 


MATRIX IS 


,//> 


, TCNT ,£46> 



1 65 . 
166 . 

167 . 

168 . 

169 . 

170 . 

171 . 

172 . 
172 . 1 
172.2 

172. 2 
1 72 . < 

172.5 

172 . 6 

173 . 

174 . 

175 . 

175 . 1 

176 . 

177 . 
176 . 

179 . 

180 . 
181 . 
182 . 
182.1 
182 . 2 

162.3 
1 82.4 
182 . 5 

162.6 
1 83 . 

184 . 

185 . 
196 . 

1 87 . 


2030 FORMAT! * l OPE N IGOP EIGENVALUES AND EIGENVECTORS........*) 

I F 1 IOL • NE « 0 ) 

1 CALL CNORM! CWR , CW I , ACL»NS ,F VEC ,VECRN, VEC IN, VECR , VECI ,VRV, 

1 VIV ,VRRV,VRIV) 

GO TO 501 
46 WRI T6 ( 6,1060 ) 

1060 FORMAT ( *OFA ILURE IN HQR2* ) 

501 IF(T.EQ.O) GO TO 16 
C 
C 

C IF THE SYSTEM IS CONTIN. PRCGRAM 'CONV* WILL 
C DISCRETIZE IT. 

C 

C 

5 CALL C0NV(BAK,GA,Q,NG,GA7,GAT1,-GL,T,NC1,NG1, 

lBA,G,SC,B,AQ,Wl Q , W3Q , CQ ,Q CC ,FTQ , G2 1 , NS , NC , NNC , 

2NCNC, AT ,LT , MT , B I , F BGC ,Y , FT , K2 1 , FT I A , FT I Y, CT , QD , ICC ) 

GO TO 18 

16 CALL MAKE (BAK, BA, NS, NS) 

18 NC-NC 1 

NG=NG1 

CALL Z0T1 < G I ,NS ,NS ) 

DO 17 1=1 , NS 

17 G I C I, I )= 1 . 

IF (NC.EQ.O) GO TO 30 
C 
C 

C PROGRAM * INNER' CALCULATES THE OPTIMAL REGULATOR, 

C PROGRAM » INNEK'CALCULATES TFE S.S. KALMAN FILTER 

C 


20 CALL INNER! NS, NC»N2»A,RM,WR,WI,X.PRG»CWP»CWI, ICNT,D, I NT, 

1ACL,SC ,GN,B , B I , GL ,FBGC, CM, BA, C, Cl ,CT,MT ,LT , AT »W21 ,TCB,W11 ,VF, 
3 Y,FT, FTIA,FTIY,FTI YA, Q21.AQ21 , ICC , 

2VEC,RVEC,VECRN,VECIN,VECP,VECI,VPV,VIV,VRPV ,VR IV ,MH^ » I A ) 

IF <NZ. EO.O) GO TO 40 



168 - 


1 88, 20 I FI T. FC.O) GO TO 33 

169. C AL L I NNEK ( N S , N S , N 2 , A , R N , W* , W I , X , P RC , C WR , CW I , I C N'T , D * I N T , 

190. lQDtGNf P tRI t GI .F BGCK. GK f E/kK f CtCI fCT«HT*LTf AT .W21.TCB .Wilt VF , 

191* 310, IQN,Y, FT, FT I A, FTI Y , F TI YA ,H , N Z , 

192. 2VEC,RVEC,VECRN,VEC IN*VECR,VECI, VKV,VIV, VRRV.VR IV,NHS, IA) 

193. GO TO 39 

194. 33 CALL I NNEK( NS ,NG * N2 , A ♦ P M t WR «W I * X , PRO , CWR ♦ CW I , ICN T, 0 , I NT , 

195. 1Q,GN,R,RI,GA,F6GCK,GM,RAK,C,CI,CT,PT,LT,AT,W21 ,TCB,WU,VF , 

196. 3IC,IQN,Y,FT,FTI A,FTIY,FTIYA,H,NZ, 

197. 2VECfRVEC v VECRNv VEC 1N» VECRtVECI tVRVfVIV, VRRV.VPIVfMHS. I A) 

19P. 39 !F( INC.EG.l) GO TO 51 

199. 40 IF (IO.FQ.n GO TC 51 

200. GO TO 60 

2 Cl. 51 CONTINUE 

2C2 . C 

2 03. C 

204. C FROM HFRf TO 60 CALCUL. CF CLOSED LOOP RESPONSE 

l 205. C WILL BE OCNE \ COV . MATRIX OF THE STATES.) 

2 06. C 

2C7. C 

208. C (X-P) = ACL* <X-P )*ACLT + ( M- P ) 

2C9, C M IS BAK » P IS FT , H SHOULD BE ZERO (LATER) 

210. C ACL IS CLOSED LOOP TRANS. MATRIX. T IS TRANS. 

211. C FOR RESP.OF CLOSED LOOP ONLY (NO PRCCESS NOISE) ....M-P=Q CP CD 

212. C 

213. C 

214. I GN=1 

215. IF (NZ.GT.O) GO TC 105 

216. NZ=NZ + 1 

217. 105 CALL ZCT1 <H , NZ , NS ) 

218. DO 106 1 = 1, NZ 

219. 1C6 R ( I , I ) = 1. 

220. IF (INC.EQ.O) GO TC 200 

221. ONE = 1 . 

222. CALL AUD(ONF. »BAK,-CNE,FT,OD»NS»NS) 

223. CALL MAKE. ( SC » FT , MS ,NS ) 




224. 



225. 



2 26. 



227. 


o 

228. 

3C0 

is 

229. 


fe 

2 30. 

200 

^ fart 

231. 


0 ^ 
O 

2 32. 


& til 

233. 


& 

234. 



235. 



2 36. 

500 


237. 



238. 



239. 



240. 

400 

i— i 

241. 

401 

O) 

242. 


? 

243. 

C 


244. 

C 


245. 

C 


246. 

C 


247. 

C 


2 48. 

C 


2 49. 



250. 



251. 

81 


2 52. 



253. 



2 54, 

82 


2 55. 



2 56. 



2 57. 

83 


2 58. 



2 50. 

SC 


CALL INNEK(;NS,NS,A2,A,RM,WP,WI,X,PR0,CWR,CWI, ICNT,0, INT, 
1QC»GN*P,RI,GI»F BGCK, GM ,ACL»C»CI»CT,MT,LT»AT,W21,TCB»W11,VF 
3IC, IQN,Y,FT,FTI A,FTIY,FTIYA,H,NZ, 

2 VEC, R VEC, VECRN, VEC IN , VECR , VEC I » VR V , VI V , VR RV ,VR IV ,MHS, I A) 
CALL ADC(GNE,ACL ,ONE,SC,ACL,NS,NS) 

GO TO 400 

I F{ T. EQ .0 ) GO TO 500 

CALL INN£MNS,NS,N2,A,RM,WR,WI,X, PRC,CWR» CWI , ICNT,D, INT, 
1QC,GN,R ,RI,GI,FBGCK,GM, ACL,C,CI,CT,MT,LT,AT,W21,TCB,W11,VF 
3IQ,IQN,Y,FT,FTI A,FTIY,FTIYA,H,NZ, 

2 VEC, R V EC, VECRN, VEC IN, VECR, VEC I ,VRV, VI V, VR RV ,VR I V »MHS, IA) 

GC TO 400 

CALL INNEKC NS,NG, N2,A ,P M , W R ,WI , X , PRO ,CWR , CW I , I CNT , 0 , 1 NT, 
iq,GN, R, RI ,GA ,FBGCK,GM,ACL»C ,C I,CT ,MT,LT ,AT , W2 1 »TCB»W11,VF, 
3I0,IQN, Y,FT,FTI A,FTIY,FTIYA ,H,NZ, 

2 VEC ,R VEC, VECRN, VEC IN, VE CR, V EC I , VR V ,V IV , VRRV »VR IV , MHS ♦ I A ) 
WRITE 16,4011 

FORMAT ( 'O', • THE CLOSED LOOP STATE COV. MATRIX IS *,//) 
CALL SP ITU ACL, NS, NS) 


E< U, UT ) = C < X-P) CT 

= C( X )CT •• . FOR EXACT MEAS. 

( C=F BGC ) 

DO 81 1*1, NS 
DO 31 J = 1 , N C 
G ( I tJ )=FBGC( J, I ) 

IF < INQ.EQ.O ) GO TO 82 

CALL ADC(ONE,ACL,-CNE,SC, ACL, NS, NS > 

CALL MULT(FBGC,ACL,Q21,NC,NS,NS) 

CALL MULT(021,G ,B,NC,NS ,NC > 

WRITE(6,83) 

FORMAT! 'G',' THE CONTROL COV. MATRIX IS ',//! 
CALL SPIT! ( B »NC ,NC- ) 

RETURN 



260 , 
261. 
2 62. 
263. 
2 64 . 

265. 

266. 

267. 

268. 

269. 

270. 

271. 

2 72. 
2 73. 
2 74. 

275. 

276. 

277. 

278. 

279. 
2 80. 
281. 
2e2. 

283. 
2 84. 

285. 

286. 
287. 
2 88 . 
289. 
2 90. 

291. 

292. 
253. 
29*. 
255. 


END 

SUBROUTINE CGNV< BAK,GA,C,NG fGAT ,GAT]*GL»T , NCI ,NC1, 

1 BA,G,SC,b,AC,WlG,W3Q, CO ,CC C ,F TQ, Q2 1, NS , NC , NNC , 

2NCNC, AT,LT,NT,BI , F 6GC , Y , FT , W2 l , FT I A, FT I Y , CT ,QD , I CC ) 
IMPLICIT P£AL*8 (A-H,0-Z ) 

D I MENS I CN FTtNS ,NS) ,W2UNS,NS >,FTIA<NS,NS 1 ,FTIY(NS,NSI, 
1 AO(NNC,NNC> , W1Q < NNC,NNC ) » W3C< NNC ♦ NNC I ,CQ( NNC ,NNC I , 

1 F TQ ( NNC , NNC ) , C2 1 ( NC,NS) , GA ( NS t NG ) , CD (NS , NS \ , 

1CTC NS ,NS> * BAKIN St NS) » GL(NS»NC> ,GATING,NS) t 
1GAT1(NG*NSI »Q(NG»NG),BA(NS»NS)»G(N!S»NCI ,SC(NS,NS), 
1AT(NCNC),LTINCNC»»MTCNCNC) , 

18(NC,NCI» QOQCNNCtNNCI , 81 <NC ,NC> ♦ FBGCiNC ,NS 1 , Y< NS , NS > 


NC= NC 1 
NG=NG 1 

CALI MAKE<eAK,6A,NS,NS) 

I F( ICC. EQ.O > GO TO 110 
C 
C 

CPRFPARE NEW MATRICES 

C 

C 

CALL Z0T1 (AQtNNCtNNC) 

DO 21 1=1, NS 
DC 21 J=1 , NS 

21 AQ( I , J) =8A< I ,J> 

DO ?2 1=1, NS 

DC 22 J = 1 » N C 

22 AQ( I , J + NS ) = G( I , J > 

CALL Z0T1 (QDQ, NNC, NNC ) 

DO 23 1=1, NS 

DO 23 J=1 , NS 

23 OCQU ,J >=SC<I,J ) 

DO 24 1=1, NC 

DO 24 J=1 ,NC 

2* Q CO ( NS+ I, N S + J )= b ( I , J ) 

CALL FATQ(T, AQ, Hi Q iW 30, CQ , N NC , CCO , FTQ > 



2 96. 


DO 25 1=1 ♦ NS 

297. 


DC 25 J=1,NS 

298. 

25 

8 A< I , J>=WlQ( I.J ) 

299. 


00 26 1=1, NS- 

800. 


DC 26 J=1,NC 

301. 

26 

G( I , J ) = QDQ( I ,NS + J > 

3 C? • 

C CALCULATE 61 = I NV. CF CDQ ( N $+ 1 ♦ NS + J ) 

303. 

C 


304, 


DO 27 J=1,NC 

305. 


DO 27 K=1,NC 

306. 


I=( J-l ) 4MC+.K 

307. 

27 

AT ( I ) =QOQ < N S*K» NS + 

3C8. 


CALL MINV(NCNC,AT,NC,CE,LT,m 

309. 


00 31 J = 1 . NC 

310. 


DO 31 K=1,NC 

311. 


I=( J-1)*NC=«-K 

312. 

31 

BI(K, J)=AT( I) 

313. 


DC 23 1=1 *NC 

314. 


DO 28 J = 1 »N S 

315. 

28 

Q21 ( I.J ) = ODG(NS + I , J ) 

316. 


CALL MULT ( B I »Q2 1 , F BGC ,NC .NC ,NS) 

317. 


CALL MULT (G.FBGC.Y .NS . NC . NS ) 

316. 


DO 29 I =1 . NS 

319. 


DO 29 J=1.NC 

3 20. 

29 

GUI, J)=W1Q( I.NS + J) 

321. 


CALL MULTt GL .FBGC . FT.AiS .NC .NS ) 

322. 


DO 30 1=1, NS 

323. 


DC 30 J=1 »NS 

3 2*- . 

30 

SCt I , J )=QDQ( I* J » 

325. 


ONE =1 . 

3 26. 


CALL ADC< ONE , SC ,-ONE , V, SC , NS.NS) 

327. 


CALL ADC ( CNE.3A .—ONE. FT .BA.NS.NS) 

328. 


IF (NG.EQ.O) GO TO 140 

329. 

110 

IF ING.EQ.O) GO TO 120 

3 30. 


DO 60 1=1 .NS 

3 31. 


00 60 J=! ,NG 



3 32. 

333. 

334. 
3 35. 
3 36. 

337. 

338. 

339. 
3 40. 

341. 

342. 

343. 

2 44. 
345. 

3 46. 
347. 
248. 
349. 
3 50. 

351. 

352. 

353. 

354. 
3 55. 
356. 
3 57. 
3 58. 
5 59. 
2 60. 
361. 
3 62. 

363. 

364. 

365. 
3 66 . 
3 67. 


G AT ( J » I ) =GA ( I , J I 

cO CONTINUE 

CALL MULTIGt GATfGATl'NGfNG*NS > 

CALL MULT(GA,G4T1 ,QD,NS ,NG,NS) 

GG TO 130 

120 CALL Z0T1 ( Q 0 » NS ,NS ) 

130 CALL EAT(ICC,T,BAK,Y,W21,FTI4,FTIY,CT,NS,Q0,FT) 

IF ING.EQ.O* GO TO 131 
WR I T F ( 6 ,7 20 > 

72C FORMAT ( 'O' , * THE CISCRETE FRCCFSS NOISE COVARIANCE MATRIX : '//) 

CALL SP!T1(QD,NS,N'S) 

131 CALL MAKEfBAKt Y,NS, NS ) 

IF (NC.EC.O GO TC 140 
I F ( ICC .EG. 1 ) GO TC 140 
CALL MA KE ( B A »Y * NS » NS > 

CALL MULT(W21,G,GL ,NS ,NS,NC ) 

WRITE (6,811 

81 FCPMAT <'0',‘ THE DISCRETE CONTROL CISTR. MATRIX IS... ',//) 

CALL SPITlIGLtNS ,NC) 

140 RETURN 
END 

SUBROUT IN E INNE K ( NS, NC, N2 , A ,RM, WR,WI ,X , PR 0 »CW R , CW I , ICNT, q 9 I NT , 

1 SC,GN,B,BI,G,FBGC,GM, EA,C,Cl,CT,MT,LT, AT, W21, TCB ,W1 J. ,VF , 

3 I G, IQN, Y, FT, FT I A, FT! Y ,FTI Y A ,H ,NZ , 

2VEC,RVEC, VECRN, VEC IN , VE CR , V EC I , VR V , V I V , VF P V ,VR I V , MHS , I A I 
IMPLICIT PEAL *8 (A-H, 0-Z ) 

DIMENSION A { N 2 , N 2 ) ,R M ( N2 , N2 I , WR (N 2 I , WI ( M2 > , X( N 2 , N2 I , PR Cl N 2, N2 ) ♦ 

1 C WR (NS) ,CVM (NS) ,ICNT( N2 I ,01 N2 I , I! T (N2 } , BA(NS,NS) , 

2GN( NS, NS) , P VEC (NS ,NS ) , VECRN( NS ,NS ) , VEC IN( NS,N S > , V£CR< NS , 

7NS ) , 

8 VEC !( NS, NS) , VRV ( N S » , V IV ( NS ) ,VRR V ( N S J , VR I V ( N S > , W1 1 ( NS ,N$ ) , C ( N S > , 
9LT(MHS) ,MT(MHS) ,B(NZ,NZ) ,GMN$»NSI» 

3W?1(N$,N$ ) ,TCB( N2 ,NS> ,AT( MHS> ,C I (NS) ,CT (NS ,NS> ,B1(NZ,NZ ) , 

5Y(NS, NS), FT (NS, NS) ,FT IA(NS,NS )» FT IVINS* NS) , FT I YA ( NS , NS ) ,H(NZ,NS) , 
4G(NS,NC) , FBGC ( N S , N Z > , VEC (NS ,NS ) , V F ( N2 , N 2 ) , SC ( NC , NC ) 

REAL*4 FMT(20) 



2 6 8 . 
369. 
37C. 

371. 

372. 
3 73. 
37A. 

3 75. 

376. 

377. 

378. 

379. 
3 80. 

381. 

382. 

383. 
386. 
3 87. 
3P8. 
3 89. 
3 SO. 
3 91. 
3 92 . 
293. 

394. 

395. 
3 96. 
3 97. 
3 98. 

3 99. 
400. 
401 . 

4 02 . 

403. 

404. 
4 05. 


C 

C***OUTPUT OPTION'S 

C IOL =1 IP THE OPEN LCCP El GENS YSTEM IS DES IR ED — TOThFRWISE ITL = 0 

IOL = 0 

C*** ORDER CF THE SYSTEM DYNAMICS 

C ( NO TE : THE DIMENSIONS CF THE SYSTEM DYNAMICS »CONTRCLt 

C MEASUREMENT , AND CCST MATRICES MUST DE SPECIFIED. 

C IN THE DIMENSION STATEMENTS ABOVE) 

M = 2* NS 

MH = NS o 

N = M 

C*4*0RDER OF THE CONTROL SYSTEM 

C NC * NUMBER OF SEPARATE CONTROLS 

NCC = NC 

C 4 4#0RD FR OF THF: MEASUREMENT SYSTEM 

C ( NOTE : I F ONLY THE CONTROL PROBLEM IS TO BE SOLVED THEN N0=C1 

NO = C 
NG= 1 


CALL Z0T1 IB I ,NZ ,NZ > 

OC 23 I ■= l * NZ 

23 BI( I, I) = 1/8(1 ,1 ) 

C***CALCULA T ICN OF KALMAN G A I NS : FORM A T I ON OF THE HAMILTONIAN 
C 


*4 

44 

* FT I 

* 

* 

* Y* FT I 

44 


4,4 
* 4 


FT I * A 


F + Y*F T I 4 a 


4 * 


C*4444*4*44*444444*44 444 44 4*: 
500 CONTINUE 
00 24 I = 


***£ U 8A NSXNS OPEN LOOP STATE 
***G IS NSXNC STATE DISTP. MTRX. 
*** e IS R NZXNZ MEASUR. NOISF MTR 
***SC IS Q NCXNC STATE NCISE MTR 
*** H IS NZXNS MEASLR..T IS TRANS. 
*4*1 1$ INV» A=HT*R I 4H Y=C*0*GT 
***THE EG. APE: X( N + l > = F*X(N ) + G«K 
ZIN)-H*X(N)+LfO=F(K4KTI , ( AVFRAGF ) 
R=E(L*LT) 

44#4*4*44444444*44444#4$4#4444:*!}c4:((]{: 


1 ,NC 




406. 


DC 24 J = 1 t MH 


407. 


PRO ( I , J ) = 0 .DO 


4Cfi. 


DG 24 K = 1 * NC 


4 C9 . 

24 

PRO ( I * J ) = PRO(I.J) + SC( 1»K)*G(J ,K> 


410. 


DC 25 I = l»l*H 


4H. 


DO 25 J = 1,MH 


412 . 


Yd t J ) =0. CO 


413. 


On 7 5 K - 1 , 5C 


414. 

25 

Y(I ,J)=Y( I , J )+&( I ,K)*PRO(K,Jl 


415. 

c 

CALCUL, A , A=H T#R I*H 1*1=0!) SCFPY TENF R C P A PY PT]Y 


416. 


00 26 1*1. NZ 


417. 


00 26 J=1 » PH 


416. 


ocfj( I t J 1 = 0. DC 


419. 


DC 26 K = i,rz 


4 20 • 

26 

PRC ( r , J )*PRQ( I . J > + BI ( I f K) *H<K , J ) 


4 ?1 . 


DO 27 I-l.Mh 


422r 


OP 2 7 J -1 » MH 

1 

423 • 


FTlYC 1, J)=O.PO 

h-» 

<1 

424, 


DC 27 K=1,NZ 

1 

* 23. 

27 

FT I Y ( I , J ) = F T I Y ( I.Jl + HCKtn^FPOlK.J) 


426. 

C 

calculate ft 


427. 


00 3001 1=1. MH 


4 2£ • 


CG 3001 J=1.MH 


4 29. 

2 C01 

F T ( I * J ) =0. DO 


4 30. 


DC 3007 1=1 »MH 


4 31. 


DO 3 002 J=l,f»H 


4 32. 

3 CC2 

FTC I , J )=6A( J , It 


4 33. 

C 

CALCUL A TE FTINV 


4 ?4 . 


DC 227 J= 1 , NS 


435. 


DC 327 K=?,NS 


4 36. 


I=( J-l )=rNS + K 


4^7. 

327 

Ann =ftc k» j i 


43 f o 


CALL MINV(MHStAT» M F.DEtLT, M) 


4 ?9 * 


DC 32 P J=lt';S 


4 4 6 « 


32 5 K=1,N$ 


441 . 


1=1 J-l )*NS+K 



r a « 




i 

i-* 

U1 

I 


<*42. 

443. 

444. 

445. 

446. 

447. 
446. 

449. 

450. 

451. 

452 . 

453. 

454. 

455. 

456. 

457. 

458. 
4 59. 

460. 

461. 
462 • 
*63. 

464. 

465. 

466. 

467. 

468. 

469. 
4 70. 

471. 

472. 

473. 

474. 
4 7*. 
416. 
477. 


328 FT(K»J)=AT<n 

C MINV R FPL ACE FT 8Y FT IN V, SO FROM HERE CN : F T=FT I NV 

C CALCUL ATE FT INV*A= FT I 1A , A = HT*R I *H 

DQ 3003 1 = 1, MH 
CO 3003 J = 1 » MH 
3C93 FTl A( I , J)=0. DO 

CALL' CM'PRUlHHSf MHS.MHS, FT , F TI Y ,FT I A ,MH , MH , MHJ 
C CALCULATE Y*FTI NV=FTI Y 

CALL GMPROtMHS, MHS ,MHS, Y ,FT ,FT I Y, MF, MH, MH ) 

C CALCULATE Y*FTIA=FTIYA 

CALL GMPRD(MHS , MHS, MHS, Y , F TI A , FT I YA ,MH , MH , MH > 

C DETERMINE THE SUBBLOCK RM22 

DQ 3004 1 = 1, MH 
DO 3004 J=1,MH 

3004 RM( I + MH»J+MH)=B A( I ,J)+FTIYAU f J) 

C CALCULAE SUBBLOCK RM21 

DO 3005 I =1 , MH 
DO 30C5 J=1 » MH 

3005 RMU + MH,J > = FTIY ( I , J) 

C CALCULATE SUBBLOCK RM12 

DC 3006 1-1 »MH 
DO 3006 J= 1 , MH 

3006 RM( I , J+MH)=FTIA( I , J) 

C CALCULATE SUBBLOCK RM11 

DO 3007 1=1, MH 
DO 3007 J=1 » MH 

3007 RM1 !, J)=+FT< I,J ) 

DO 3008 J=1 , M 

DO 3006 1=1, M 

3008 A(1 , J )=RM( I , J) 

$***$*$*****#** *** ***** *********** 

CALL BALANC < M, N , A ,LOW, I H IGF, D ) 

CALL ELMHES i M, N , L CW , I H I GF , A , I NT ) 

CALL FQR2 ( M ,N, LOW , I H IGF, A , WR , W I , X , ICNT , £46 ) 

CALL ELMBAK \ M, LOW , I H IGH, N, A , I NT , X ) 

CALL BAL8AK ( M, N » LOW , IH IGH ,N,D , X » 



47b. 

479. 

480. 
461. 

482. 

483. 
4 9*- . 
485. 
4 86 . 

487. 

488. 

489. 
4S0. 
491. 
4 92. 
4 93. 

4 94. 

495. 

496. 

497. 

498. 

499. 

5 OC . 
5C1. 
502, 
5 03. 
5 04. 

505. 

506. 

507. 

508. 
5 09. 

510. 

511. 

512. 

513. 


C************ ##*##*** St******************** 

CALL G4INI N ,NS ,NC,RM,WR,WI ,X,GN»fcll , T CB, 

1W21 ,4T,LT,C,CI,CT, PHS,MT» 

CALL MAKF<6A,GN,NS,NS ) 

IF ( IQN.EQ.l ) GO TO 3e9 
WRI TE ( 6 ,1 401 ) 

1401 FCRMATC '0* i ’THE MATRIX M (THE CCV. OF EST. ERROR BEFORE MFASUR.) 
1 IS •,//) 

CALL R A P R N T ( MH, MH , W , 8, GN,4 , • ( 3 (1 X , 1 PD 1 3. 6 ) )• ) 

C CALCULATION OF P , P = F INV ( M-Y I FTI NV , FI N V= FT I NVT 
00 650 1=1, MH 
DC. 650 J= 1 , MH 

650 GNU, J)=GN< I ,J)-Y( I, J> 

OO 651 1=1, MH 

DC 651 J=1 , MH 
PRO ( I , J)=C.OC 
DO 651 K=1,NH 

6 51 PRC ( I , J >= PRO < I , J ) + GN { I,K)*FTtK,J) 

DC 652 1=1, MH 
DO 652 J=1,MH 
GN ( I » J ) =0.00 
DO 652 K= 1 , MH 

6 52 GM I, J)=GN( I , J ) +FT (K , I) *PRC (K , J ) 

W R I T E ( 6 ,1402 I 

CALL MAKE ( FT »GN ,NS , NS ) 

1402 FORMAT ( * THE MATRIX P ( THE CCV. OF EST. ERROR AFTER MEASUR.) 

1 IS ',//) 

CALL R A PR NT (MH,MH»MH,8» GN, 4 , •( 8 ( 1 X , 1PD1 3. 6 » )» ) 

C 

C CALCULATION OF FILTER GAIN 
C 

C CALCUL. FILTFR GAIN K=GN*HT*RI 

4 FORMAT 1 8( 2X , 1PD 13.6 I) 

DO 800 I = l ,MH 
00 800 J = 1 , NZ 
PRO! I , J) =0.00 



177- 


514. 

515. 

516. 

517. 

518. 

519. 

520. 

521. 

522. 

523. 
5 24. 

525. 

526. 

527. 

528. 
5 29. 
5 30. 

531. 

532. 

533. 

534. 

535. 

536. 

537. 

538. 
5 39. 

540. 

541. 

542. 

543. 

544. 

545. 

5 46. 

547. 

548. 

549. 


CO 800 K = 1,NZ 

800 PRO(I,J) = PRO(T,J) + HK,1)*SI(K,J) 

C F6GC=KALMAN=K=GN*PR0 

DC 801 I = I,ks 
DC 801 J = 1,NZ 
FBGCilfJ) = 0.00 
DO 801 K = 1,NS 

801 FBGCUtJ) - FBGC ( I » J ) + GN( I ,K ) *PRC (K, J) 

WRITE! 6,977) 

577 FORMAT ( ' ' , 'THE FILTER GAINS K = P*HT*RINV ARE:*,//) 
DC 968 I « 1 ,NS 

568 WRITE (6,4) (FBGC(I,J),J = 1,NZ) 

WRITE ( 6 ,5 998 ) 


CALL RGAI N(M,NS,NC,RM,V)R,V(I,X,GN,Wll ,TCB, 
1W21 ,AT,LT ,C,CI #CT , NHS » MT ) 


5998 FORMAT ( • 1 * , • ESTIMATOR EIGENVALUES AND E I GENV ECTCRS. . . • ) 

1VIV,VRRV^VR1V» M { C,CI?CT .ns.rvec.vecrn.vecin.vecr.vecuvrv 

1400 CONTINUE 


» 


GO TO 389 

46 WRITE(6,1060) 

1060 FORMAT ( *0 FA I LJR E IN HCR2M 
339 RETURN 
END 


1 !r| BR Sr T JEV2y E i t i!;? C,h2,A,RM ’ WR,M » X » PRO * CWR » CWI *ICN T ,0 T INT, 

ACL,SC,GN,B,BI,G,F BGC , GM,BA,C,CI,CT,MT,LT,AT,W21,TCB,Wli,VF, 
3Y,FT,FTIA,FTIY, FT I YA , G2 1 , AQ?1 , ICC, 

^IMPLICIT* REAL*8 I A^ ^ r ^ ^ I * VR V *VI V , VR RV , VR I V »MH S, I A > 

, SlJJfUf ! CN A tN2,N2),RM(N2,N2),WR(N'2),WI(N2),X(N2tN2),PROCN2,N2)t 
} : CWHNS) * ICNT< N2 >,0 (N2) , INTIN2) , BA ( NS, NS ) , ACL ( NS, N S) , 
2GNMS,NS> ,RVEC(N5,NS> ,VECRN (NS , VS ) , VEC I N< N S ,NS ) ,VECR(NS, 

7NS) ,021 (NC,NS), AQ21(NC,NS) , 

8 VEC It NS, NS) ,VRV(NS) ,VIV(NS) ,VRRV(NS) ,VRIV(NS> , W11(NS,NS ) , C (NS >, 
LT ( MHS ) ,MT (MHS ) ,B( NC , NC ), GM (NS ,NS J , 

3W21 (NS,NS ) ,TCB( N2 ,NS) ,AT(MHS) ,C I(NS) »CT (NS.NS ) .81 (NC.NO . 



I 


? 


550. 

551. 

552. 
553 . 

554. 

555. 

556. 

557. 

558. 
5 59. 

560. 

561. 

562. 

563. 

564. 

565. 
5 66. 
5 71. 
5 72. 
5 73. 

574. 

575. 

576. 

577. 


5Y(NS, NS), FT (NS, NS ) t FT IA < NS, NS ) , FT l YCNS , NS ) , FT I YA ( NS ,NS > , 

4G< N$,MC ) ,FBGC<NC ,N$) , VEC ( N S ,NS ) , VF (N2 , N2 > , S C t NS , NS > 

REAL*4 FMT ( 20 ) 

C 

C#**0UTPUT OPTIONS 

C I0L=1 IF THE OPEN LOOP EICENSYSTEM IS DESIRE C — TOTFER to IS E IOL=0 

I CL = 0 

C *** ORDER OF THE SYSTEM DYNAMICS 

C ( NO TE i Tl-E DIMENSIONS CF THE SYSTEM DYNA MI CS .CONTROL , 

C MEASUREMENT .AND COST MATRICES MUST BE SPECIFIED 

C IN THE DIMENSION STATEMENTS ABOVE) 

Ml = 2*NS 
MH * NS 
N = M! 

c ***orq ER QP THE CONTROL SYSTEM 

C NC = NUMBER OF SEPARATE CCNTROLS 

NCC = NC 
NO = 0 
NG= 1 

C ********** ****** ****** ****** ****** ***** 

C IN I TAL I Z AT I ON OF MATRICES 
IF ( ICC.EO.i ) GO TO 500 
DO 7 1= 1 » NC 
DO 7 J = 1 » NC 



-6L T 


r 


578, 


7 0I( I , J| = 

0. DO 

579. 


DO 23 I = 

1 » NC 

5 80. 


23 8111,11 = 

1 ✓ B C I , I) 

sex. 

582. 

C***CALCULATION 

C 

OF CONTROL GAINS 

5 83. 

C 

** 

AA 

58A. 

C 

* 

A 

5 85. 

C 

* F+Y*FTI*A 

-Y*FTI A 

5 86. 

C 

* 

A 

587. 

C 

* 

A 

588 . 

C 

* 

A 

5 85. 

c 

♦ -FT I * A 

FT! A 


FORMATION OF CONTROL HAMILTONIAN 


*** F ANO FT ARE THE OPEN LOOP 
DYNAMICS MATRIX AND TRANSPOSE 
IS NCXNC CONTROL WEIGHTING 
MAT RI X 

IS THE NSXNS STATE WEIGHTING 
MATRI X 

*** FT I = FT INVERSE f Y=GM*BI*GM 
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5 90. 
591. 

5 92. 

593. 

594. 

595. 

596. 

597. 

598. 

599. 

600 . 
601. 
602 . 

603. 

604. 

605. 

6 06 . 

I 607. 

608 . 

609. 

610. 
611 . 
612 . 
6 13 . 

614. 

615. 

616. 

617. 

618. 
619. 
6 20 . 
621. 
622. 
623. 
6 2 *. 
6 25. 


C ** ** ***GM IS THE NSXNC CCNTRCL 

C OISTP I8UTI0N MATRIX 

C * ** #***♦#*#****!*#*** « ***$***#**♦#*$*********$+#*********:#**** 

500 CONTINUE 

DO 24 I = 1 , INC 
DC 24 J = 1 ,MH 
PRO (I »J) = 0.00 
DO 24 K = 1 * NC 

24 PROII,J» - PRO(I,J) + B I ( I ,KI*G(J ,K» 

00 25 I - l,MH 

DO 25 J = l »MH 
Y' I » J )=0. DO 
DO 25 K = l » NC 

25 Y f I ,J»=Y< I,J)4GU ,K)*PRC<K f J) 

C CALCULATE FT 

DO 3001 I =1 ,MH 
CO 3001 J=l,MH 

3001 F T ( I , J) =0.00 
DO 3002 1 = 1, MH 
DC 3002 J = l » MH 

3002 FTI I , J>=BA< J,I> 

C CALCULATE FTINV 

DO 32 7 J=1,NS 
DO 327 K = 1 1 N S 
I=<J-1>*NS+K 

327 ATI I >=FT( K, J) 

CALL MINV (MHS,AT,MF,D£t LT,MT) 

DO 328 J= 1 , N S 
DO 328 K= 1 , NS 
I = I J— 1 * *NS+.K 

328 FT(K,J)=ATI I) 

C MINV REPLACE FT BY FTINV, SC FROM. HERE ON : FT = FTINV 

C CALCULATR FT INV *A = FT I I A ,A= SC 

DC 3003 1=1, MH 
00 300 3 J = 1 » MH 

3003 FTI All , J»=O.DO 



reproducibility of the 
ORIGINAL PAGE IS POOR 


626. 

627. 

626. 

629. 

630. 
631. 
622. 
6 33. 

634. 

635. 

636. 
627. 

638. 
629. 
640. 
6 41 . 

642. 

643. 

644. 

645. 

646. 

647. 

648. 

649. 
6 50. 
651. 
6 52. 

653. 

654. 

655. 
6 56. 

657. 

658. 

659. 

660. 
661 . 


CALL GMPRDl MHS, MHS ,MHS, FT , SC, FT I A , PH, MH , MH > 

C CALCULATE Y*FT I NV = FT I Y 

CALL GMPRDl MHS, MHS ,MH$, Y,FT *FT I Y » Ml-, MH, MH ) 

C CALCULATE Y*FT I A= FT I YA 

CALL GMPRD< MHS , MHS, MHS ,Y , FTI A , FT I YA , MH , MH , PH ) 

C DETERMINE THE SUBBLOCK RM11 

DO 3004 I =1 , MH 
DO 3004 J = 1 1 MH 

3 004 RM( l , J ) =-BA ( I ,J )-FTIYA( l,JJ 
C CALCULAE SUBBLOCK PM12 

DO 3005 1=1, MH 
DC 3005 J = 1 »MH 

3005 RM( I , J+MH ) = FTIY ( I , J) 

C CALCULATE SUBBLGCK RM21 

DO 3006 1 = 1, MH 
DC 3006 J=1 * MH 

3006 RM( I+MH,J|=FTIA(I,J) 

C CALCULATE SUBBLGCK RM22 

DO 3007 1=1, MH 
DO 3007 J=1 , MH 

3007 RM( I + MH,J+MH)=-FT( I, J ) 

DC 3008 J=1 , M 

DO 3008 1=1, M 

3008 AU ,J )=-RMC I ,J| 

C ******************** ********************* 

CALL BALANC ( M, N» A ,10 W, I HIGH , D » 

CALL ELMHES C M, N , LOW , IH IGH, A, I NT I 

CALL HQR2 (M,N,LQW»IHIGH,A,WR,WI,X,ICNT,&46> 

CALL ELMBAK ( M, LOW , I H I GH, N , A, I NT, X I 

CALL BALBAK C M, N , LOW , IH IGH , N, D , X ) 

C *********** ********* *********** ********** 

CALL RGAI MM,NS»NC»RM,'dR,wl,X,GN»W!l,TCB, 
1W21,AT,LT,C,CI,CT , MHS ,MT) 

WRITE! 6,1 401 ) 

1401 FQRMAT ( *0 * , • THE RICATTI GAIN MATRIX*,//) 

CALL RAPRNT ( MH, MH , MH , 5, GN , 4 , • ( 5( IX, 1P01 3. 6) ) • ) 



662 . 

663 * 

664 . 

665 . 

666 . 
6 67 . 
668 . 

669 . 

670 . 

671 . 

672 . 

673 . 
6 74 . 

675 . 

676 . 

677 . 
6 7 3 . 

t m 679 . 

“ 6ec. 

681 . 
682 . 

683 . 

684 . 

685 . 

686 . 
687 . 
6 98 . 

689 . 

690 . 
6 51 . 

692 . 

693 . 

654 . 

655 . 

656 . 

657 . 


C CALCULATION OF FEEDBACK GAIN 
C 

C CALC. FEEDBACK GAIN U=-BINV*GT7FTINV(GN-A) , A=SC 

4 FGRMAT(8(2X, 1PD13.6H 

C CALCULATE GT 

DO 81C 1=1, MH 
DO 810 J=i,MH 

810 PRO(I»J> = GN( I,J)-SC< I »J> 

DO 811 I = 1,MH 
DO 811 J = It MH 
GN( If J ) =0.00 
OG 811 K = 1,MH 

811 GN( I . J I = GN(ItJ) + FT C I ,K l + PROC K ,JI 

00 300 I = 1,NC 

DO 800 J = 1 f MH 

PkOi If J ) = 0.00 
DC 600 K = 1 , MH 

800 PROUtJ) = PRQdtOI ♦ G( K , I } *GN ( K, J ) 

DO 801 I = 1 » NC 

DO 801 J = 1 ,NS 
F BGC ( I , J) = O.DO 
DO 801 K = 1 ,NC 

801 FBGCC I » J) = FBGC(ItJ) +B I( I ,K ) *PRO (K , J ) 
c ***thE OPTIMUM FEEDBACK CONTROL CAINS 

CO 802 I = 1 »NC 

DO 802 J = ItNS 

802 FBGC(IfJ) = -FBGCUtJ ) 

IF (ICC.EQ.l) GC TC 145 
WRITE(6,977> 

577 FORMAT ( • *,«THE CONTROL CAINS ARE:*, //I 
00 96 8 I = IfNC 

568 WP ITE ( 6 ,4 ) (FBGCUfJ)tJ = 1,N$> 

C***THE CLOSED LOOP DYNAMICS MATRIX 
145 CONTINUE 

DO 150 I = 

DO 150 J = 


1 * NS 
1 ,NS 



658. 

6 <39. 
700. 
70] • 

702. 

703. 

704. 

705. 

706. 

707. 
708* 
709. 

7 10. 

711. 

712. 

713. 

714. 

715. 
71b. 

717. 

718. 

719. 

720. 

721. 

722. 

723. 

724. 

725. 

726. 

727. 

728. 

729. 
7 30. 

731. 

732. 

7 2 3. 


PRO! I , J > = 0. DO 
DO 150 K = 1,NC 

150 PRO!I,J) = PRO( I ,J KG! I ,K)*FBGC (K , J) 

DO 160 I = 1,NS 
DC 160 J = 1 , NS 

160 ACL< I t J)=BA! I,J >+PRO( I, J) 

IF ( ICC.fc'C.O ) GO TO 803 
C C=CPRI ME-6I *Q21 

CALL MULT! 81 ,Q2 1 , AG21 ,NC,NC ,N$ ) 

ONE =1 . 

CALL ADQ( ONE »FBGC ,-ONE, AQ21 ,FBGC,NC,NS) 

WRITE(6»957I 

587 FORMAT!* ♦ , • THE CONTROL GAINS ARE:*,//) 

DO 988 I = 1 ,NC 

538 WRI TE < 6 ,4 ) ( F 8GC <1 , J ) ♦ J = l,NS) 

803 WRITE! 6, 170) 

170 FORMAT ( *0 ' , * THE CLOSED LCCP DYNAMICS MATRIX IS..*,//) 

DO 180 I = 1 , NS 

180 WRITE (6,4) (ACL(I,J),J = 1,NS) 

WRITE! 6,5998) 

5558 FORMAT! *1*, 'CLOSED LOCP EIGENVALUES AND EIGENVECTORS...*) 
1VIV V^R VR*i V)^ ^ CtCI.CT »NS»RVEC»VECRN»VECIN»VECR,VECI»VRV» 

1400 CCNTINUE 
GO TO 389 
46 WRITE(6,1060) 

1060 FORMAT! ’CFAILURE IN HQR 2 * ) 

389 RETURN 

END 


C 

c 

c 

c 

c 

c 


SUBROUTINE E ATQ (T,A,Wl,W3,C,N,UC,FT) 


THIS SUBROUTINE COMPUTES THE TRANSITION MATRIX AND IT*S INTEGRAL. 
THE StRIES IS TRUNCATED WHEN THE LARGEST ELEMENT OF THE LAST TERN 
THE SERIES IS LESS THAN l.E-03 TIMES THE SMALLEST ELEMENT CF THE 
SERIES. WRITTEN BY R. MAINE 8/17/71 



734. 

735. 
7 36. 
737. 
7 38. 

739. 

740. 
7 41. 
7 42 . 

743. 

744. 

745. 

746. 
7*7. 
74 P • 
7 49. 
750. 
7 5 i • 

752. 

753. 
75*. 

755. 

756. 

757. 


C 

C COMPUTES PHI < T/ 3)=EXP<A*T/P ) AND ITS INTEGRAL 

C USES HOLLOWING RELATIONS 3 TIMES TC OBTAIN PHI(T) AND OS INTEGRAL 

C PHI (?T » = PHIC TM<*2 

C INTEGRAL ( PHI (?.T )) = INTFGR AL ( PHI ( T))*( I + PHI( Til 

C 

C INTEC-R AL ( PHI ( 2T ) «2T 1 = 1 NTEGRAKPHI ( T )«T |* ( i+PEI ( T ) + 

C + I NT EG P AL ( PH I IT ) l*PHI< T I *T 

IMPLICIT REALMS! A— h,0-Z ) 

DIMENSION A(N,N ) , M. ( N ,N ) , W3 l N » N ) *C (N*NI » 

1QC( Nf N> »FTl N,N> 

CO 5 I = I » N 
DC 5 J=1,N 
5 FT( I, J »=A( J, I> 

Z=T/8. 

CALL MULTfODtAfC ,M,N,N) 

CALL MULT (FT.C. W1 ,N,N,N) 

CALL MULTCFT.QD ,W3 ,N,N, M 
Z2=Z*Z/2. 

CALL ACG(ZtQD,Z2. »C»QD,N,N) 

ON E = 1 . 

CALL ADD< ONE,QD ,Z2 ,W3,QC, N, M 
Z3=Z2*Z*?./3. 

CALL ACC(GNE,QD ,Z3 ,W1 


( “C. 

759. 

760. 

761. 

762. 

763. 

764. 

765. 

766. 

767. 

768. 
7 69. 


NX = N 
N T = 2 4 

CALL Z0T1 ( W1 »N, N ) 

DO 1 1 = 1, N 
WI( I » I ) =1 . o 
CALL NAKF (W3,W1,N,N) 
G =1 . 0 

W3MAX = 1 .E+50 
T=T/g. 

00 7 1=1, NT 
BB= I 

G=G*T/BB 



185 


7 70. 

771. 

772. 

773. 
7 74. 

775. 

776. 

777. 

778. 

779. 

7eo, 

781. 

782. 

783. 
7 64. 

785. 

786. 

787. 

788. 

789. 

790. 

791. 

792. 

793. 

794. 
755. 

796. 

797. 

798. 

799. 

800. 
801. 
802. 
803. 

8 04. 
805. 


3C 


40 


WlMIN*l.E+50 
DO 30 J = 1 f N X 
DO 30 K=1 ,NX 
IF< Wl( J.K). NE.O. 

CONTINUE 


•AND. CABS (W1 < J , K ) ) .1 7. W1 M I N ) W 1M I N=DABS < W 1 ( J , K|) 


W3MAX1 =W3MAX*T/ BB 
CALL MULT ( A , W3 ,C ,N» N , N ) 

CALL MAKE ( W3 ,C . N, N) 

W3MAX=0. 

DC 40 J-l .NX 
DO 40 K=1,NX 

IF CDABS( W3( J ,K U *GT . W3MAX) |W3MAX=0ABS( W3(J,K) ) 
CCNTI NLE 


W3MAX=W3MAX*G 

CALL ACC (QNE,W1,C,W3,W1,N,N) 

IF< W3MAX .LT. W1MIN*1. £-03 »GC TO 70 
7 CONTINUE 

WRITE (6.1000) W1MIN.W3MAX.W3MAX1 
1 000 FORMAT! 'OERROR IN EAT * , 5X , ■ WlM I N = NE15.6, 5X, • N3MAX 
117X,*W3MAX1 = ' » E15.6 ) 

70 CONTINUE 
DC 90 K=1 ,3 


= * »E15.6/« 


DO 71 L=1,N 
DC 71 J=1,N 

71 FT( L» J)~W1( J.L) 

CALL MULT (Q0.W1 ,W3 .N. N, N) 

CALL MULT(FT,W3 .C.N.N.M 
CALL ADD< OME.QD , ON E , C ,Q0 , N , N) 

CALL MA KE ( W3.W 1, N ,N) 

CALL MULT LW 1, W 3 , C ,N, N, N ) 

CALL MAKE (Wl.C.N.N) 

T=T*2. 

90 CONTINUE 

WRITE (6,51) I 

51 FORMAT ( '0 *, *THE EXTENDED TRANSITION MATRIX*, 15,' TERMS*/) 
CALL SPIT! (Wl.NX.NX) 


1 


t 



3C9. 

BIO. 

811 . 

812. C 

313. C 

8 14. C 

815. C 

8 16. C 

817. C. 

619. C 

319. C 

8 20. C 


RETURN 

END 

SUBROUTINE EAT ( I C C , T t A , Wl , W2 , W3 , W* ,C , N ,Q0 , FT ) 

THIS SUBROUTINE COMPUTES THE TRANSITION MATRIX AND IT’S INTEGRAL. 
THE SERIES IS TRUNCATED WHEN THE LARGEST ELEMENT OF THE LAST TERM 
THE SERIES IS LESS THAN l.E-03 TIMES THE SMALLEST ELEMENT OF THE S 
SEF. IPS. 


COMPUTES PHI (T/8> = EXP(A*T/3 ) AND ITS INTEGRAL 

USES FOLLOWING RELATIONS 3 TIMES TO OBTAIN FHKT) AND ITS INTEGRAL 


6 21 . 
822. 
623. 
9 24. 
325. 
6 26 . 

327. 

328. 
6 29. 
8 30. 
3 31. 
8 32. 
833. 
8 34, 

835. 

836. 

837. 
8 38. 

839. 

840. 

841. 

842. 

843. 
8 44. 


C PHM2T)*PHl(T)**2 

C INTEGRAL (PHI ( 2T >> = INTEGRAL ( PHI (7)>*< I + PWKTM, 

C 

C INTEGR AL< PHI (2T> *?T) = INTEGRAL! PHUT )*T>*< I +PHI (T> + 

C *INTEGRAL(PHI(TI )*PH I ( T ) *T 

IMPLICIT REAL*8 I A-W,0-Z I 

DIMENSION A ( N f N » » Wl( N ,N » , W2 <N ,N > , W3( N ,N ) , W4 (N , M ,C < N, N » , 
1 QD I N » N ) »FT ( N » N) 

DO 5 1 = 1, N 
DC 5 J = 1 , N 
5 Wl( I , J )=A( J , I ) 

Z=T/8. 

CALL MULT( QD,W1 , W2,N,N,N) 

CALL MULTI A ,W2, W1 ,N,N,N) 

CALL MLLT(A,QO, W3 ,N,N,M 
Z2=Z*Z/2. 

CALL ACC(Z,QD,Z2,W2tQC,N:,NI 
ONE =1 . 

CALL ADC( ONE * CD * Z2 ,W3,0C,N,M 
Z3=Z2*Z*2./3. 

CALL ADCt ONE,QD , Z3,W1 ,QC,N,M 
NX =N 
N T= 2 *■ 

CALL ZCT1(W1,N,N) 



— £81 


845 . 

846. 

847. 

848. 

849. 
8 50. 

851. 

852. 
8 53. 
8 54. 
855. 
8 56. 
8 57. 

858. 

859. 

860. 
861. 
862. 

863. 

864 . 
e65. 
866 . 

867 . 

868 . 

869. 

870. 

871. 

872. 

873. 
8 74. 
8 75. 
8 76. 
877, 
8 78. 
8 79. 
880. 


CALL MAKE (W2,Wl,N,N> 

DC 1 1 = 1, N 
1 Wl< I, I >=1.0 

CALL MAKE (W3,W1,N,NI 
G = 1 .0 

W3MAX=l.E+50 
T=T/8 • 

00 7 1=1, NT 
8B= I 

G=G*T/B8 
WlM IN=i .E + 50 
W2MIN=l.E+50 
DO 30 J = 1 , N X 
DC 30 K=1 , NX 

I F( Wl ( J ,K ) . NE.O • .AND. CAB S C Wl ( J , K > ) .L T. W1 M I N > WlMI N=D ABS ( W1 ( J ,K > > 

1 F ( W2 (-J ,K ). NE.O . .AND. CABS(W2(J,K)) .L T . W2MI N J W2M I N=DAB S ( W2< J , K ») 
30 CONTINUE 

W3MAX1=W3MAX*T/ BB 

CALL ADD (ONE»W2,G,W3»W2,N,M 

CALL MULT { A,W3 ,C ,N,N,N) 

CALL MAKE CW3,C,N,M 
W3M AX=0. 

DC 40 J=1 , NX 
DO 40 K=1,NX 

IF COABSC W3( J ,K ) > .GT. W3MAX ) W3M AX=DAB S { K3 < J , K ) ) 

40 CONTINUE 

W3M AX =W 3MAX * G 

CALL ADD (0NE,Wl,G,W3,Wl,N,M 

IF(w3MAXl .LT. W2M IN*1. E-03 .AND. W3MAX .LT. WlM I N*l. E -03 > GO TC 70 
7 CONTINUE 

WRITE ( 3,1000) WlMIN,Vi3MAX,W2MIN,W2MAXl 
1000 FORMAT! ’OERROR IN EAT • , 5X , « Wl M I N. = » ,E1 5 . 6 , 5 X, • W3M AX =',S15.6/» •, 
117X,'W2MIN = * ,E 15. 6,5X, 'W3MAX1 =*,E15.6> 

70 CONTINUE 
DO 90 K=1 , 3 
DO 71 L = 1 , N 
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B81. 
882, 
8 83, 
884* 

885. 

886 . 

887. 

888 . 

889. 

890. 

891. 

8 92. 

893. 

894. 

895. 

896. 
8S7. 

898. 

899. 

900. 

901. 

902. 

903. 

904. 

905. 

906. 

907. 

908. 

909. 

910. 

911. 

9 12. 
913. 
9 14. 
915. 
9 16. 


DC 71 J = 1 r N 
71 FT(L. J)=W1( J.L) 

CALI MULT(GD,FT,W3,N,N. N) 

CALL MULTIHl»W3.FT.NtN,M 
CALL ADC( ONE ,QD .ONt.FT, CC.N.N) 

CALL MAKE (W3,WI,M,N> 

CALL MULT ( m 1 , W 3 . C ,N , N , N ) 

CALL MAKE (Wl.C.N.M 
CALL MUL T (W2.W3»A»N»N»N) 

DC 80 J=1,N 
80 W3< J, J >=W3( J.J) +1. 

CALL MULT (W2.W3 ,C , M.N.N) 

CALL MAKE <W2,C,N,N> 

T=T*2. 

90 CONTINUE 

IF ( ICC.EO.l ) GO TO 91 
WRITE (6.51) I 

51 FORMAT ( *l * , * THE TRANSITION MATRIX*. 15,* TERMS'/) 
CALL SPIT1 (Wl.NX.NX) 

91 RETURN 

END . 

SUBROUTINE MULT( A.B.C.N.M.K) 

IMPLICIT REAL*8 ( A-H.O-Z ) 

DIMENSION AIN. Ml .C(N.K) 

DC 10 1=1, N 
DO 10 L=1,K 
XX=0. DO 
DC 11 J = 1 » M 
11 XX = XX + A( I ,J )*B( J.L ) 

10 C (I ,L » =XX 
RETURN ‘ 

END 

SUBROUTINE SPITHA.N.M) 

IMPLICIT REAL*8 ( A-H.O-Z ) 

DIMENSION A ( N. M I 
CO 10 N 1= 1 * N 



BEPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


917. 

918. 

919. 

920. 

921. 

922. 

923. 

924. 

925. 

926. 

927. 

928. 

929. 

930. 

931. 

932. 

933. 

934. 

935. 

936. 

937. 
93fl. 
9 39. 
9 40. 

941. 

942. 

943. 

944. 

945. 

946. 

947. 

948. 

949. 

950. 

951. 

952. 


10 WRITE(6,20) < A< NI , N J I ,N J=1 , M 

20 FORMAT (8( IX, 1PD13.6) ) 

RETURN 

END 

SUBROUTINE Z0T1<A,N,M) 

IMPLICIT REAL*8 I A-h,0-Z ) 

dimension a { N,M ) 

DO 10 I =1 , N 
DO 10 J=1,M 
10 A ( I , J)=G.DO 
RETURN 
END 

SUBROUTINE MAKE ( A , B, N , M I 
IMPLICIT REAL*8(A-H,0-Z ) 

D I MENS I CN A ( N,M ) ,e(N,N) 

DO 10 1=1, N 
DO 10 J =1 , M 
10 A ( I , J ) =BC I, J) 

RETURN 

END 

SUBROUTINE ADD ( X , A, Y , B , C , N ,M ) 

IMPLICIT REAL*8 ( A— H» Q-Z ) 

DIMENSION A(N,M ) ,B(N,M) ,C( N,M) 

DO 10 I =1 , N 
00 10 J=1,M 

10 C (I ,J) = X*A<I ,J) +Y*BC'I ,J ) 

RETURN 

END 

SUBROUTINE R APR NT ( NMA X , M , N , L, A , IDI M ,FMT ) 
REAl*8 AINMAX,N» 

DIMENSION FMTIIDIMJ 
NU= L 

DO 20 NL=1,N,L 
IF (NU.GT.N) NU=N 
DO 10 I = 1 , M 

10 WRITEI6,FMT)(.A{ I,J),J = NL,NU> 



I 


<D 

? 


953. 
9 54. 

955. 

956. 

957. 

958. 

959. 

960. 

961. 
96?. 
963. 
9 64. 

965. 

966. 

967. 

968. 

969. 

970. 

971. 

972. 

973. 
9 74. 
9 75. 

976. 

977. 

978. 

979. 
9 60. 

981. 

982. 
993. 
9 64. 

985. 

986. 

987. 
9 88. 


WPITE(6,100 ) 

100 FORMAT! • * ) 

20 NU=NU+L 

RETURN 

END 

SUBROUTINE CDIV < A ,B, Ct C, E, F ) 

IMPLICIT RE AL*8 <A-H,C-Z> 

C 

C THIS SUBROUTINE COMPUTES THE COMPLEX DIVISICN 

C 

c t + f*i * t a ♦ b*i )/<c ♦ c*n 

C 

T-C*C*D*D 
E = ( A*C+9*0> /T 
F=( 9*C-A*DI/T 
C 

RETURN 

END 

SUBROUTINE HQR2 ( MM, N,L OW ,H I ,H, HP , V! , F , ICNT ,* ) 

C 

IMPLICIT REAL*8 (A-H,C-Z) 

DIMENSION H(NM,NM| f WR(NMI»WI(NM),F(NM,NMJ,ICNT(NM» 
INTEGER HI ♦ HU . EN» EL 
REALMS IM f IV t MA ( 50 hMB( 501 
DATA EPSM / Z 34 100 CO 0 00 00 CO 00 / 

LOGICAL LAST 
C 

DO 101 1=1, N 
DC 9101 J=I,N 
F (I ,J)=O.ODO 
9101 F(J,I)=O.ODO 

101 F( I .1 ) -=1.000 
C 

IF ( N.GT.2)MA(3 J=H<3, 1) 

IF (N.LT.4) GO TO 5102 
CO 102 I = 4 » N 



989, 



MA< I 1 =H< I ,1-21 

990. 


102 

MBi I > = Ht I, 1-3 1 

991. 

C 



992. 


9102 

10WM1=L0W-1 

993. 



HIl=HI +1 

994. 

C 



995. 

C 


STORE KNOWN EIGENVALUES 

996. 

C 



9 97. 



IF U0WV1.LT. 1) GO TO 9103 

998. 



DC 103 1=1, LOWM 1 

999. 



WR { I ) =H i I , I > 

1000. 



WI( 11=0 .000 

l 

1 001 . 


" 10 3” 

Torn i > =o 

1002. 

C 



1003. 


9103 

IF (HU.GT.N) GO TC 9104 

1004. 



DC 104 I=HI 1 , N 

1005. 



WRi I ) = Hi I *-I > 

1006. 



WIi I ! = 0.0D0 

1007, 


104 

ICNTi 11=0 

1008. 

C 



1009. 


9104 

E N=HI 

1010. 



T=0.0D0 

1011. 

C 



1012. 

C 


DETERMINE EIGENVALUE EN 

1013, 

c 



1014. 


105 

IF <EN.LT.LCW1 GO TO TOO 

1015. 

c 



1016. 



I TS =0 

1017. 



NA=EN-1 

1018. 

c 



1019. 

c 


SEARCH FCR SPLIT 

1020. 

c 



1021. 


200 

L = EN 

1022. 


1C6 

IF < L.EQ.LOW > GC TO 109 

1023. 



IF (DABS (H(LfL-l) ) • LF.EPSM’MDAESi H(L-1»L-1 1 1 

1024. 


X 

+ DABS< H( L ,L 1 ) ) ) GO TC 109 



1025, 


L = L-1 

1026. 


GO TO 106 

1027. 

C 


1028. 

C 

TEST FOR CONVERGENCE 

1C29. 

C 


10 30. 

109 

X = H ( EN » EN ) 

1031. 


IF (L.EQ.EN) GC TO 110 

1 032. 


Y=h( f\ A tNA ) 

1033. 


< EN » N A ) *H 1 N A » EN > 

1C34. 


IF IL.F0.NAI GO TO 111 

1035. 


IF I ITS. EQ. 301 PETURN1 

1036. 

C 


1037. 

C 

COMPUTE SHIFT 

10 38. 

C 


1039. 


IF ( ITS. NE. 10. AND. ITS.NE.PO) GC 

1040. 


T =T + X 

1041. 


DO 112 I=L0W,EN 

1042 • 

112 

H(I ,II=HI I ,1 i-x 

1043. 


S=D A BS i H ( E N* NA) >+ D AE S < H ( NA , EN- 2 

i 0 44 • 


Y=. 75DO*S 

1045. 


X = Y 

10 46. 


W= -0 .437 5D04 S* S 

10 47. 

C 


1048. 

c 

OP ROTATION 

1 G49. 

c 


10 50. 

113 

ITSMTS+l 

10 51. 


EL=NA-L 

1052. 


DO 114 MM=l,Fl 

1053. 



1054. 


Z=H( 

10 55. 


R -X-£ 

10 56. 


S = Y-Z 

10 57. 


P=(R«S-W)/H( R+ltM)+H(H t M+l> 

1058. 


C=H(M+1,M+1 I-Z-R-S 

1059. 


R=H( M + 2 » M +1 ) 

1060. 


S = OABS( P> +DABS (CI+CABSt R ) 



-061 


1061. 


P=P/S 


1062. 


Q=Q/S 


1063. 


R=R /S 


1064. 


IF (M.EQ.L) GO TO 115 


1065. 


IF 1 DABS( H ( M ,M-1) 1* ( CABS < C MOABS 1 P ) ) . L E. EPSM*DABS< P ) 

1 066 . 

X 

*(OABS(H< ^-l,M-m + DAES(ZKDA8S(H(M*l t M+Ul)) 

1067. 

X 

GO TO 115 


1068. 

114 

CONTINUE 


1069. 

C 



1070. 

115 

MP 2 = M+ 2 


1071. 


00 116 I = MP2,EN 


1072. 

116 

H( I » 1-21=0.000 


1073. 


MP 3=M+ 3 


1074. 


IF (FN.LT.MP3) GO TO 9117 


1075. 


DO 117 I=MP3,EN 

* 

1076. 

117 

H ( I » 1-31=0.000 


1077. 

9117 

CONTINUE 


1078. 

C 



1079. 


DO 118 K=M »N A 


1080. 


LAST=K.EQ.NA 


1081. 


I F ( K.EQ, Ml GO TO 119 


1082. 


P=H( K,K-1 ) 


1083. 


C=H ( K+ 1 »K -1 ) 


10 84. 


R =0.000 


1085. 


IF ( .N0T.LAST»R=HIK+2,K-1» 


1086. 


X=DABS< P) +0A 8S ( C) +DAeS(R ) 


1 0 87 . 


IF (X.EQ.O.OCO) GO TO 118 


1088. 


P = P/X 


1089. 


Q=Q/X 


1090. 


R=R/X 


1091. 

119 

S=O.ODO 


1092. 


IF (P.EQ.O.OCO) GO TO 1119 


1093. 


IF (DLGGIOI CABS (PI) .GE.-38.ODOI 

S=P*P 

1094. 

1119 

IF (Q.EQ.O.OCO) GO TO 2119 


1095. 


IF <DLQG1Q<CABS<Q )) .GE.-38.GO0) 

S=S+Q*Q 

1096. 

2119 

IF (R.EQ.O.OCO) GO TO 3119 




194- 


1097. 


IF IDL0G10ICABS IR 1 ) .GE.-38.OCO) S=S+R*R 

1098. 

3119 

S=DSQ9T(S ) 

1099. 


IF IP.LT.0)S=-S 

1100. 


IF ( K.EQ. GC TO 120 

1101. 


HIK.K-1 »=-$*X 

1102. 


GC TO 9121 

1103. 

120 

IF (L.NE.M>H<K t K-l>=-HK,K-l) 

1104. 

9121 

P=P + S 

1105. 


X=P/$ 

llOo. 


Y = G/S 

1107. 


Z=P/S 

1108. 


G=G/P 

1109. 


P =R / P 

1110. 


DC 121 J= K t N 

11 11. 


P=HIK, J)+G*HIK+1,J J 

1112. 


IF (LAST) GO TC 122 

1113. 


P=P+R#H{ K + 2, J) 

11 14. 


H(K+2»J» = H(K+2tJI-P*Z 

1115. 

122 

H< K + l, J)=H(K+1 ,J) -P*Y 

1116. 

121 

HI K,JI =H 1 K, J l-F*X 

1117. 


J -MI MOIENt K+3) 

1118. 


DO 123 1*1, J 

1119. 


P» X*HI I ,K )+Y *H ( I ,K4l > 

1120. 


IF ILAST ) GO TO 124 

1121. 


P=P+Z*H( I t K+2) 

1122. 


HI I ,K*21 = F{ I r K42)-F*R 

1123. 

124 

Hi. I »K+U = HI I ,K + 1 I-P*Q 

1124. 

123 

HI I . K ) =H I I t K )- F 

112 5. 


DC 125 I = LOW .HI 

1126. 


P=X*F( I r K 1+ Y*F I I ,K+1 » 

1127. 


IF ( LAST1 GO TC 126 

1128. 


P=P + Z*FI I tK+2) 

1129. 


FI I ,K*2 I- F I I »K+2)-F*R 

1130. 

126 

FI I ,K+ l)=Ft I t K + 1 )- F*Q 

1131. 

125 

FI I ,K> =F< I » K )— P 

1132. 

lie 

CONTINUE 



1133, 

C 


1134. 

C 

END OF QR ROTATICN 

1135. 

c 


1136. 


GO TO 200 

1137. 

c 


1138. 

c 

CNE REAL ROOT IS DETERMINED 

11 39. 

c 


1140. 

110 

H( EN ,EN l = X+T 

1141, 


WR!EM=H<EN,EN> 

1142. 


W I { EN) =0.000 

1143. 


IF ( EN.NE.l IH(EN.NA) =0.0C0 

1144. 


I CNT { EN ) = ITS 

1145. 


EN=N A 

1146. 


GO TO 105 

1147. 

c 


1148. 

c 

- — TWO ROOTS ARE DETERMINED 

1149. 

c 


1150. 

111 

P=( Y-XJ/2.0D0 

1151. 


Q=P*P+W 

1152. 


? = DSCRT( DABS <01 ) 

11 53. 


X=X + T 

1154. 


H(FN,EN)=X 

1155. 


H{ NA »MA) =Y+T 

11 56. 


IF {NA.NE.1)H(NA,NA-U = 0.0D0 

1157. 


ICNTIEN) =-ITS 

11 58. 


ICNT (NA) = ITS 

1159. 

c 


1160. 


IF K.LE.O.OOO) GO TO 2C1 

1161. 

c 


1162. 

c 

ThO REAL ROOTS 

1163. 

c 


11 64. 


IF < P. LT. 0.000) Z=-Z 

1165. 


2= P + Z 

1166. 


WR( NA ) =X + l 

1167. 


s=x-w/z 

1168. 


WR CENI =S 



-961 


1 It**, 




Wi INM=U.UD0 

1170. 




Mil EN)=O.ODO 

1171. 




X=H< ENtNA) 

1172. 




r=dscrt< X*X+Z*Z ) 

1173. 




P-X/P 

1174. 




Q-l/P 

1175. 




DO 203 J=NA,N 

1176. 




Z =H ( NA • J) 

1177. 




HtNA.J >=Q*Z*P*H(EN, J) 

117b. 


20 3 


H(ENt J»=0*H4 EN, J)~P*Z 

11 79. 




CO 204 1=1, EN 

1160. 




Z=H <1 ,NA> 

1181 . 




H ( I ,NA>=Q*Z+P*H{ I , EN ) 

11 62. 


204 


H ( I ,EN)=Q*H< I,ENI-P*Z 

1183. 




CO 205 l = LGW ,H I 

1164. 




Z=F ( I » NA ) 

11 85. 




F ( I ,NA)=Q*Z+P*F( I , EN ) 

1186. 


205 


F<I,FN)=Q*F< I,EM-P*Z 

1187. 




H ( FN ,NA)=O.ODO 

11 88. 




GO TC 202 

1189. 

C 




119C. 

C 




1151. 

C 



TWO COMPLEX RCCTS 

1192. 

C 




1193. 


201 


WP l NA ) = X + P 

1154. 




WR(EN)=X+P 

1195. 




WI(NA)=Z 

11 56. 




Wl» EN) =-Z 

11 57. 

c 




1158. 


2C2 


EN =E N-2 

1199. 




GO TC 105 

1200. 

c 



EN C OF EICENVALUF ITERATION 

12 01. 


100 

RNR M = C. CD 0 

12 02. 



K* 

1 

1203. 



DC 

210 1=1 , N 

12 04. 



DO 

209 J=K , N 



1 2 C5. 209 RNRM=RNRM+DABS( H( l , J ) ) 

1206. 210 K= I 

1207. C 

1208. C DETERMINE THE EIGENVECTORS OF THE TRIANGULAR MATRIX STORED 

1209. C IN H ANC OVERWRITE THEM ON H 

1210. C 

1211. EN=N 

1212. 400 IF (EN.LT.2) GO TO 401 

1213. C 

1214. P=WR(EM 

1215. Q=WI CEN) 

1216. NA = EN- 1 

1217. C 

1218. IF (G.NE.O.ODOl GO TO 21 2 

1219. C 

1220. C REAL EIGENVECTOR CORRESPONDING TO REAL EIGENVALUE 

1221. C 

i 1222. I =N A 

w 1223. 206 IF U.LT.U C-0 TO 220 

V 1224. C 

1225. 5 = H(ItENI 

1226. IP1 = I+1 

1227. IF (NA.LT.IP1) GO TO 9214 

1228. DC 214 J= IP1 »NA 

1229. v. 214 S=S+H( I r J >*H(J t ENI 

1230. 9214 Y =0 .000 

1221. IF II.NE* 1)Y=HI I, I-i) 

1232. 215 Z=P-H(I,n 

1233. IF (Z • EQ. O.OD0 ) Z=E PSM<RNRM 

1234. IF ( Y.NE. O.OGO) GC TG 9213 

1235. C 

1236. H ( I »EN ) =S /Z 

1237. GC TO 213 

1238. C 

1239. 921? 1=1-1 

1240. P=HM,EN) 



1241, 


I P2= 1+2 

1242. 


IF (NA.LT.IP2) GO TC 9216 

1243* 


OQ 216 J=IP2,NA 

1244. 

216 

R-R+HII.J >*H<J»ENJ 

1245. 

921 fc 

1 1 ,n-P 

1246. 


X-H(I,I+1 ) 

1247. 


IF I DABS! W ) .IE. DA 0S f Y 1 > GO TC 217 

1246. 

C 


1249. 


PK=Y/ W 

1250. 


Z=< S-RM*R 1 /< Z+RM*X) 

1251. 


H ( I + 1 f ENI =Z 

1252. 


H( I ,EN >=( -R-X+Z »/R 

12 53. 


GC TO 213 

12 54. 

C 


12 55. 

217 

RW=W/Y 

1256. 


XM RM*S— R )/( X + R M* Z ) 

1257. 


H{ I +1 i EN ) =X 

12 58. 


HU,5N) = (-StX*ZI/Y 

12 59. 

C 


1260. 

C 


1261. 

213 

1 = 1-1 

1262. 


GO TC 206 

1263. 

C 


126*-. 

220 

H( FN,FN}=1.000 

1265. 


GO TO 211 

1266. 

C 


12 67. 

C 

COMPLEX E IGF NVECTCR CORRESPONDING TO COMPLEX EIGENVALUE 

1 2 6P. 

c 


12c9* 

212 

I = N A 

12 70. 

299 

IF (I.LT.l) GO TO 35C 

1271. 

C 


12 72. 


E =H ( I » EN ) 

1273. 


S = 0. 0 DO 

12 74-. 


I PI = 1 + 1 

12 75 . 


IE < IP1.GT.NAI GO TO 9301 

12 76. 


ne 3ci j = i p l »na 



“66 1 


1277. 


R =R+H ( I » J ) *H ( J t NA ) 

1278. 

301 

S=S+H( I ,J >*H( J,EM 

1279. 

9301 

CONTINUE 

1280. 


Y=O.ODO 

1281. 


IF ( I ,NE» 1 ) Y=HI I,I-1» 

1282. 

302 

z=hii , n-p 

1283. 

C 


1284. 


IF (Y.NE.O.OCOI GC TG 9303 

1285. 

C 


1286. 


CALL CDIV(-P»-Sf2'QfH(I f NA) f H(I,EN) ) 

1287. 


GC TO 303 

1288. 

C 


12 89. 

9 303 

I =1-1 

1290. 


R A=H ( I «EN> 

1291 . 


SA=O.ODO 

12 92. 


I P2=I + 2 

1293. 


IF (NA.LT.IP2) GO TC 5304 

12 54. 


DO 3 04 J= I P 2 , N A 

1295. 


R A= RA +H( I ,JMHI J,NA 1 

1296. 

304 

SA= SA +H( I r J ) *HI J, EM 

12 97. 

9 304 

CONTINUE 

1298. 


V*=HM , I >-P 

1299. 


X = H ( 1,1+1) 

13 00. 


IF (DA8S( W ) + DA6SI O.LE. DABSIVI ) GO TO 305 

1301. . 

C 


1302. 


CALL CDIV (Y »C.OCO,WtG»RK, INI 

1303. 


R=R— RM*RA + I^*SA 

1 304. 


S=$-RM*$A-IM*RA 

1305. 


Tl=RI**X-2 

1306. 


T2=IM*X-Q 

13C7. 


CALL CDlVIRtS*Tl>T2 v RV»IV) 

1308. 


7i=-R A-X* R V 

13 09. 


T2=-$ A-X* I V 

1 5 1 0 . 


CALL CDIVITl ,T 2 , W ,Q»H(I , N A ) ,HI ( ,EN) ) 

1311. 


GO TO 306 

13 12. 

C 




200 - 


13 13. 


'J . 


G V. = w/ Y 

1314. 




I M = 0/ V 

1 3 \ r . 




p A = r.A -94!* R+ I M* S 

1316. 




S A= $ A-RM# S- [F*R 

1317* 




T i = o M*Z-I M* Q-X 

1 3 1 8 . 




T2= 1M*Z+RM*C 

1 319. 




CALL CDIVIRA ,SA,T1»T2,4V,IV) 

132C. 




HI ! ♦ NA ) = ( IV*Q,-RV*Z-F ) / Y 

13 21. 




HI I . c M =~ ( S+ IV* 7+FV*Q 1/Y 

132?. 

C 




1323. 


306 


H t r +1 ,N',J =R V 

13 24. 




H ( I + 1 , C N> =!V 

1325. 

C 




1 376. 

c 




1327. 


303 


1 = 1-1. 

1 323. 




GO TC 29.9 

1329. 

c 




13 30. 


3 50 


H ( ,-N,NA» =1.0 00 

1331. 




H ( r N » n N > = 0. 0 DO 

13 32. 




f N=KA 

1333. 

c 




1334. 


211 

EN= 

= EN-1 

13 35. 



GC 

T 0 409 

1336. 

c 




13 37. 

c 


— 

•— END EIGENVECTORS OF TRIANGULAR MATRIX 

1338. 

c 



’ 

13 39, 

c 




2 3 40. 


4-ai 

I F 

(w m I.EQ.O. 000 1 H 1 1 , II = 1.3C0 

1341, 



IF 

(LDWM1.LT. U GO TO 94C2 

1342. 



DC 

40? 1 = 1 r LC;*'l 

1343. 




IP1=I+1 

1344. 




UO 402 J=IP1 ,N 

1345. 


40 2 


FI I » JI=HU tJ I 

1346. 

c 




13 47. 


540? 

IF 

(HU.GT.M GO TC 404 

1348. 



DO 

40 3 I = HI 1 * f " 



201 


1349. 



IF (I.EQ.N) GO TO 94C3 

1350. 



IP1=I+1 

13 51. 



DO 403 J=IP1 f N 

1352. 

403 


F(I t J)=H(I»J> 

1353. 

C 



13 54. 

9403 

I F 

(LOW.GT.HH GC TO 404 

1355. 


00 

416 J= H 1 1 

1356. 


DC 

416 I=LOW,HI 

1357. 



2=0. ODO 

1358. 



00 405 K=LOM » H I 

1359. 

405 


Z=Z+F( I, K)*H(K, J) 

1360. 

416 


F ( I , J ) =Z 

13 61. 

C 



1362. 

404 

J=H I 

1363. 

9404 

IF 

< J. LT. LOW > GO TO 413 

1364. 

C 



1365. 



IF ( WKJI.EQ. 0.000 > GC 

13 66. 

C 



1367. 



IP=J-1 

1368. 



DO 4C8 I =LOW ,HT 

1369. 



Y=0 .ODO 

1370. 



Z=0.000 

1371. 



DO 409 K= LOW » J 

1372. 



Y=Y+F{ I » K ) *H 1 K » IP I 

1373. 

409 


Z=Z+F< I,K >*H (K» J) 

1374. 



F(I,IP)=Y 

1375. 

408 


Ft I , J)=Z 

1376. 



J= I P 

1377. 



GO TO 406 

13 78. 

C 



13 79. 

407 


00 410 I=L0W,HI 

1380. 



Z =0*0D0 

13 81. 



DO 411 K= LOW » J 

1382. 

411 


Z=7.+F( I »K J) 

13 83. 

410 


F ( I » J ) =Z 

1384. 

C 





i .1 H S . 

1386. 

1367. 

1388. 

1309. 

1390. 

1391. 

1392. 

1393. 

1394. 

1395. 

1396. 
1 397. 
1398. 
1 399. 

1400. 

1401. 

1402. 

1403. 
14 04. 

1405. 

1406. 

1407. 

1408. 

1409. 

1410. 

1411. 
14-12. 
1413. 
1 4 !<► . 

1415. 

1416. 

1417. 
l4l8o 
1419. 
14 20 . 


406 J=J-1 

GO TO 9404 
C 

C PNO EIGENVECTOR DETERMINATION 

C 

413 IF <N.GT.2)H<3 t H = NA<3) 

IF <N.LT.4>RETURN 
DO 415 1 = 4 1 N 
H ( I , 1-2 > = M A ( I > 

415 HI,I-3MB(I) 

C 

RETURN 

END 

SUBROUTINE BALANC (NM t M , A , LOW , H I , D ) 

C 

IMPLICIT REALMS <A-H t C-Zt 
DIMENSION AINM,NV ) ,f)(NM) 

INTEGER HI 
. LOGICAL NCC GNV 

DATA B . B2 / 16.0D0 t 256. OCO / 

C 

L = 1 
K S N 

c 

C SEARCH FOR ROWS ISOLATING AN EIGENVALUE AND PUSh THEM DOWN 

C 

ICO J = K 

101 IE (J.LT.l) GO TO no 

R = O.OD? 

on 102 I = 1»K 

102 R = R+D ABS ( A ( J y I ) ) 

R = R-OABSI AI J, J) ) 

IF (R.NE.O.OCO) GO TO 103 
0IK» = J 

IF IJ.EC.K) GO TO 20? 

00 201 1 = 1 »K 



I 

to 

o 

v 


1421* 



F = A( I 

1422. 



All , J) 

1423. 


201 

A( I, K ) 

1424. 



OC 202 I = 

1425. 



F = A ( J 

1426. 



At J,I| 

1427. 


202 

A ( K ? I ) 

1428. 


203 

K = K — 1 

1429. 



GO TO 100 

1430. 


103 

J = J-l 

1431. 



GO TO 101 

1432. 

C 



1433. 

C 


SEARCH FOR 

1434. 

C 



1435. 


110 

J = L 

1436. 


109 

IF (J.GT.K 

1437. 



C = O.OCO 

1438. 



DC 112 I = 

1439. 


112 

C * C+DABSI 

1440. 



C = C-DABSI 

1441. 



IF tC.NE.O. 

1442. 



D C L 1 = J 

1443. 



IF (J.EQ.L1 

1444. 



OC 211 I = 

1445. 



F = A( I , 

1446. 



At 1, J) = 

1447. 


211 

At! ,L> * 

1448. 



DO 212 I = 

1449. 



F » A( J t 

1450. 



At J . I ) * 

1451. 


212 

At L » I ) = 

1452. 


213 

L = L + l 

1453. 



GO TO 110 

14 54. 


111 

J = J + l 

1455 • 



GO TO 10° 

1456. 

C 




A{ K T I) 
F 


SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE AND PUSH THEM 


L,K 


GO TO 213 

1»K 


A( L * I ) 
F 


LEFT 



204 - 


1457. 

1458. 

1459. 
146G • 

1461. 

1462. 

1463. 

1464. 

1465. 

1466. 

1467. 

1468. 

1469. 

1470. 

1471. 

1472. 

1473. 

1474. 

1475. 
1 4. 76 • 
1477. 
14 78. 

1479. 

1480. 

1481. 

1482. 

1483. 

1484. 

1485. 

1486. 
14 87<s 

1488. 

1489. 
14 90. 

1491. 

1492. 


113 LOW = L 
Hi ~ K 

IF (L.GT.K) RETURN 
DO 300 I = L . K 
3CC OH) = 1.000 
C 

c NOW BALANCF THE MATRIX IM ROWS L THROUGH K 

C 

302 NOCONV = .FALSE. 

OC 301 I = L,K 


C = 0. ODD 
R = O.ODO 
00 303 J = L » K 
IF (J.EQ.I) GO TO 303 
C = C+CA8S ( A ( J » I) > 

R = R+OABSI A ( 1, JM 

303 CONTINUE 
C 

F = 1. CD 0 
S = C+R 
C 

G = R/B 

304 IF (C.GE.G) GO TO 305 

F = F*B 
C = C.*B2 
GO TO 3C4 
C 

305 G = R* B 

306 IF (C.LT.G) C-0 TO 30 7 

F = F/B 

C = C/B2 
GC TC 206 
C 

307 IF < (C+RI/F. GE. 0.9506*$) GO TO 301 



205- 


14*53. 

1494. 

1495* 

1496. 

1497. 

1498. 

1499. 

1500. 

1501. 
15 02. 
1503. 
15 04. 
1505. 
15C6. 

1507. 

1508. 

1509. 

1510. 

1511. 

1512. 

1513. 

1514. 

1515. 

1516. 

1517. 

1518. 

1519. 
15 20 . 

1521. 

1522. 

1523. 

1524. 

1525. 

1526. 

1527. 
15 28. 


31-1 

312 

C 

301 

C 


C 


C 


102 

101 

c 

107 


104 


G = 1.000/F 
o< n = om*F 
NOCGNV = .TRUE. 

DO 311 J = L,N 
A< I , J» = A(I ,J)*G 
00 312 J = 1 t K 
A< J. I> = A< J 'I ) *F 

CONTINUE 


IF (NOCONV) GO TO 302 
RETURN 


END 

SUBROUT IN E BAL8 AK ( NM, N f LQW»HI ,H,D , 2) 

IMPLICIT REAL*8 (A-H, C-Z) 

INTEGER HI 

DIMENSION D(NM) ,Z(NM,NM) 

IF (LOW.GT.HII GO TO 107 
DC 101 I=LOW,HI 
S=D(I) 

DO 102 1 » M 

ZC I. J)=Z( 1,J >*s 
CONTINUE 

IF (LOW.LE.U GO TO 108 
L0WM1=LCW-1 
00 103 1 1 = 1 » LOW M 1 
I = L C k— I I 

k=o( n 

IF (K.EQ.n GO TO 103 
DO 1 C4 J — 1 1 M 
S=Z ( I » J ) 

Z(I tJ) = z< K f J> 

Z(K , J > = $ 



206- 


1C 3 


CONTINUE 


15 29. 

15 30* 

1531. 

1532. 

1533. 
15 34. 

1535. 

1536. 
15 37. 
15 38. 
15 39. 

1540. 

1541. 

1542. 

1543. 
154*. 

1545. 

1546. 

1547. 

1548. 

1549. 

1550. 
15 51. 

1552. 

1553. 

1554. 

1555. 

1556. 

1557. 

1558. 
15 59. 

1560. 

1561. 

1562. 
1 563. 
156*. 


108 IF tHI ,GE .N > RETURN 
Ih=Hl + 1 
DC 105 1= 1H » N 
K = D< I) 

IF (K.EO.I) GO TO 105 
DO 106 J=1,M 
S = Z € I . J I 
Z f I , J)=ZI K,J ) 

10 6 Z ( K » J ) “ S 

105 CONTINUE 

RETURN 

END 

SUBROUTINE ELMH ES < NM f N , K » L » A , INT) 

IMPLICIT R E AL*8 <A-H»C-Z> 

DIMENSION INTINMI T AINM,NM) 

LM1=L-1 

KP1=K+1 

IF <LM1.LT.KP1I RETURN 

DO 101 MMKP 1 » LM 1 
I = M 

X=O.ODO 
00 102 J = M,L 

IF (DABSC A< J ,M-1) J.LE.DABSt XI » GO TO 102 

X=A< J,M-1 > 

I=J 

102 CONTINUE 

INT ( M) = I 

IF (I.EQ.M) GO TO 103 
DG 104 J=MM1 ,N 



1*5 65. 

1566. 

1567. 
1566. 
15 69. 

1570. 

1571. 
15 72. 
1573. 
15 74. 
1575. 
15 76. 

1577. 

1578. 

1579. 
1560. 
1581. 
15 82. 
1583. 
15 64. 

1565. < 

1566. 
1587. 

1 5 88/ 

1589. 

1 590. 

1591. 

1592. 

1593. 

15 94. 

1595. 

1596. 

1597. 

1598. 

1 599. 
1600. 


C 


c 


c 


Y=A ( 'I t J > 

A( I , J)=A< M, J ) 

104 A ( M » J )=Y 
DO 105 J=1,L 

Y=A(J,I) 

A ( J * I )=AC J,M» 

105 A(J,MI=Y 


103 IF (X.EQ.O) 
MP1=M+1 


GO TO 


\ DO 


106 I=MP1 t L 
Y=AM,M-1 > 
IF (Y.EQ.O) 
Y=Y/X 

A I I » M-l )= Y 


GO 


101 


TO 106 


DO 107 J=M,N 

107 A < I ,J)=AI I f J >-Y*A(M, J ) 

DO 108 1,L 

108 A( J T M)sA< J,K)+Y*AfJ,I > 

106 CONTINUE 

101 CONTINUE 


RETURN 

END 

SUBROUTINE ELMBAK <NM,K,L t R t A t INT»Z) 

IMPLICIT REALMS (A-H,0-Z> 

DIMENSION A I NM »NM ) ,1 NTC NM I , Z( NM t NM ) 
INTEGER R 
LMK -L-K-l 

IF (LVK.LT.l)RETURN 
DO 101 MM=1 T LMK 
M= L-MM 
MP 1 = V+ 1 

DO 102 I=MP1,L 
X = A ( I t M-l ) 
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1601. IF (X.EQ.OI GO TO 102 

1602. DC 103 J=l,R 

1603. 103 Z(I . J » =Z < I »J )+X*Z<M, J ) 

1604. 102 CONTINUE 

1605. I = INT(M) 

1606. IF (I.EQ.M) GO TO 101 

1607. DO 104 J = 1,R 

1608. X*ZH,JI 

1609. Z < I * J )-Zl M, J ) 

1610. 104 ZIM,J»=X 

1611. 101 CONTINUE 

1612. RETURN 

1613. END 

1614. SUBROUTINE RGAI N< M , NS , NC, RH , WR , W I , VF, GN , W 1 1 ,Tf B , 

1615. 1W21,AT,LT,C,CI,CT,PHS,MT1 

1616. IMPLICIT REALMS <A-H,0-Z) 

1617. DIMENSION WR<M) »WI(M» ,VF(M, Ml ,GN<NS,NSJ ,RM<M,M) 

1618. DIMENSION Wilt NS, NS > ,TCBtM»NS) , W2HNS, NS) ,LTIMHS) »MT ( MHS ) 

1619. DIMENSION AT(MHS) 

1620. DIMENSION C ( NS > ,C I (NS > , CT<N S,NS ) 

1621. MH = NS 

1622. K = 1 

1623. KP = 1 

1624. KN = 1 

1625. 10 IF(K.GT.M) GO TO 200 

1626. *Z=WR (K)**2+WI( K)**2-l. 

1627. IF(WZ) 100,5000 ,50 

1628. 50 IF(MKKI) 80,75,80 

1629. C EIGENVECTOR FOR REAL E I GEN V ALUE , PO S ITI VE 

1630. 75 CONTINUE 

1631. KP = KP+1 

1632. K=K +1 

1633. GO TO 10 

1634. C EIGENVECTOR FOR COMPLEX El GENVALUE , PCS I TI VE REAL. PART 

1635. 80 CONTINUE 

1636. KP = KP+2 
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1673, 


CO 320 1=1 » MH 

16 74# 


00 320 J=1,PH 

1675# 

320 

W21 1 I * J ) = TCBd+MHvJ) 

1676. 

322 

CONTINUE 

1677. 



1678. 

C INVERT tfll 

1679. 


DO 327 J=1 . NS 

1680. 


00 327 K= 1 » NS 

1681. 


I=IJ-1I*NS +K 

1682. 

327 

ATI n = Wll<K,J ) 

1683. 


CALL MINVIMHS, AT , PH , CE TC » LT , PT ) 

1684. 


DO 328 J= 1 » NS 

1685. 


DC 328 K=1 , NS 

1686. 


I = ( J-l > #N S +K 

1687. 

328 

Wll I K, J )=AT ( I ) 

1688. 

C 


16 89. 

C CALCULATE THE PGAI N MATRIX 

1690. 


00 325 I L=1 , PH 

1691 • 


DC 325 JL=1,MH 

1692. 


GN( IL, JL) * 0.00 

1653. 


00 325 KL =1 t PH 

16 94. 

32 5 

GNI IL, JL)=GNIIL , JL I+W21 II L , KL ) *W1 1 IKL , J L ) 

1655. 


GO TO 6000 

1696* 

5C0C 

WRITEI6,1) 

1697. 

1 

FORMAT! »1 '»• ZER G E IGENV ALUE-R I CAT T I EQUATIGN 

16 98. 

6000 

RETURN 

1699. 


ENO 

17 CO. 


SUBROUTINE MINV I NR , A , N , 0,1 , P> 

1701. 


IMPLICIT REAL*8 IA-H,0-Z) 

17C2. 


DIMENSION AINM) »L(NM),M!NP) 

17C3. 


DOUBLE PRECISION A , 0, B I CA , HGLO 

17C4. 


0=1.000 

1705. 


NK=-N 

1706. 


DO 80 K=1,N 

1707. 


NK=NK + N 

1708. 


L(K ) = K 



ITS- 


17 G9. 



nIK »=k 

1710. 



KK=NK+K 

1711. 



B IGA=A ( KK ) 

1712* 



DO 20 J=K,N 

1713. 



I Z-N* ( J-l ) 

1714* 



DO 20 I =K *N 

1715. 



I J= IZ+I 

1716. 


10 

IF( DAeS(BIGA)- DABS ( A( I J ) ) ) 15*20.20 

1717. 


15 

B IGA= A ( IJ 1 

1718. 



L 1 K I = I 

1719. 



M (K )= J 

17 20. 


20 

CONTINUE 

1721. 

C 



1722. 

C 


INTERCHANGE RCWS 

1723. 

C 



1724. 



J=L( K ) 

1725. 



I F < J- K » 35,35,25 

1726. 


25 

K I = K - N 

1727. 



DO 30 1=1 ,N 

1728. 



K I=K I +N 

1729. 



HCLD=— A <KI 1 

1730. 



JI=KI-K+J 

1731. 



4(KI)=A(JI) 

1732. 


30 

A ( J I ) =HOLD 

1733. 

C 



1734. 

C 


INTERCHANGE COLUMNS 

1735. 

C 



1736. 


35 

I =M CK J 

1737. 



IE< I-K) 45,45,38 

1738. 


38 

JP=N« { 1-1 | 

1739. 



DO 40 J=i,N 

1740. 



JK=NK+J 

1741. 



JI=JP+J 

1742. 



HOLC=-A(JK) 

1743. 



A 1 JK ) = A 1 J I ) 

1744. 


4G 

A ( J I ) = HOLD 



I /<o. 

L 



1746. 

c 


DIVIDE COLUMN BY MINUS PIVOT (VALUE CF PIVOT ELEMENT IS 

1747. 

c 


CONTAINED IN BIGA) 

1748. 

c 



1749. 


45 

I F < B I G A > 48,46, 48 

1750. 


46 

D=O.OD0 

1751. 



RETURN 

1752. 


48 

DC 55 1=1, N 

1753. 



IF(I-K) 50,55,50 

17 54. 


50 

I K=NK+ I 

1755. 



A(IK)=A( I K ) /(-B IGA) 

1756. 


55 

CONTINUE 

17 57. 

c 



1758. 

c 


REDUCE MATRIX 

1 759. 

c 


■ 

1760. 



00 65 1=1, N 

17 61. 



I K = NK+ I 

1762. 



HOL D=A ( IK) 

1763. 



I J= I-N 

1764. 



DO 65 J=1 ,N 

1765. 



IJ= IJ+N 

1766. 



IF(I-K) 60,65,60 

1767. 


60 

IF(J-K) 62,65,6 2 

1768. 


62 

KJ= I J-I+K 

1769. 



A ( I J ) =HCLD* A < KJ I+AIIJ) 

1770. 


65 

CONTINUE 

1771. 

c 



1772. 

c 


DIVIDE ROW BY PIVOT 

1773. 

c 



1774. 



KJ=K-N 

1775. 



DC 75 J = 1 * N 

1776. 



KJ=KJ+N 

17 77. 



IF(J-K) 70,75,70 

1778. 


70 

A(KJ) =A(KJ) /BIGA 

1779. 


75 

CONTINUE 

1780. 

c 







1781. 

C 


PROOUCT OF P IVOTS 


1782. 

C 




1783. 



0=0*8 I GA 


1784. 

C 




1785. 

C 


REPLACE PIVOT BY RECIPROCAL 


1786. 

C 




1787. 



AlKKl = C1.0D0)/BIGA 


1788. 


80 

CONTINUE 


1789. 

c 




1790. 

c 


FINAL ROW AND COLUMN INTERCHANGE 


1791. 

c 




1792. 



K=N 


1793. 


100 

K = ( K- 1 ) 


1794. 



IFIK> 1 50 ,1 50 »l 05 


1795. 


105 

I=L ( K ) 


17 96. 



IFII-K) 120,120,108 


1797. 


108 

JQ=N*( K-l) 

1 

W 

1798. 



J R= N* (I -1 ) 

M 

CO 

1799. 



DO 110 J=1,N 

1 

1800. 



JK= JQ+ J 


1801. 



HOL D=A ( JK ) 


1802. 



J 1= JR + J 


1803. 



A ( J K I =— A < J I ) 


1804. 


110 

A ( J I > = HOLD 


1805. 


120 

J=M(K» 


1806. 



I F ( J- K ) 100 , 100 ,12 5 


1807. 


125 

K I = K-N 


1808. 



DO 130 1=1, N 


1809. 



K I=KI+N 


1810. 



HOL D = A< KI » 


1811. 



J I = K I -K+J 


1812. 



A(KIJ=-A(jn 


1813. 


130 

A { J I ) =HOLD 


1814. 



GG TO 100 


1815. 


150 

RETURN 


1816. 



END 
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1B17. SUBROUTINE CNORMI W Z , W Y , VEC » N$ , P VE C , VECR N , V E C I N , V FCR, VEC I , VR V, 

1818, IV IV ,VRRV, VR IV) 

18194 IMPLICIT R E A L *8 <A-H,0-Z) 

1820. DOUBLE PRECISION hZ,WY,VEC 

1821. DIMENSION WZ ( NS ) , W Y < NS > , VEC I NS , NS ) 

1822. DIMENSION RV FC< NS , NS ) T VECRN (N S »N$ ) ,VECIN< NS ,NS > , VEC R ( NS , N S ) , 

1823. 1 VECI( NS, NS) ,V RV < N S > , V I V < NS ) ,VR RV < NS ) , VR IV ( NS ) 

18 24. 7 FOR MATl • » ) 

1825. 8 FORMAT < • OREAL EIGENVALUE ('♦12,' I ',14X,*REAL EIGENVECTOR 

16 26. 1( ' , 12, ' >••♦•••• .'//IX,' l * ,F15.Q,* >+J(' ,F15.8, • )« ,13X,» < ' ,F15.8,*> « 

1827. 2) 

ie28. 9 FORMAT! 44 XfM*, FI 5 . 8, ' ) •) 

1829. 10 FORMAT < 'CCOMPLEX E IGENV ALUE ( ' , I 2 , * ) •, 14X , ’COMPLEX E I GEN VEC 

18 30. 1T0R( • ,12, * ) ...'//1X,» C ’,F15.8,’ ) + J< • ,F15.8, • ,13X, •< * ,F15.8 

1831. 2, *)+J( * ,F15.8,» >* ) 

1832. 11 FORMAT (44X, ' ( ' , F15.8, * MJ( * ,F15.8, •) •) 

18 33. C 

1834. LM=0 

1835. LP=0 

1836. LC=0 

1837. DO 999 K =1,NS 

1838. LK -K+l 

18 39. REMCD = O.DG 

1840. IF (DABSC W<Y( K)> .LT.l.D-10) GO TO 998 


1841 . 

1842 . 

1843 . 

1844 . 

1845 . 
1 846 . 

1847 . 

1848 . 

1849 . 
18 50 . 

1851 . 

1852 . 


IF CLM. EQ.l . AND. K. EG. NS I GG TO 1000 

IF (LM.EQ.l) GO TO 992 

LC=LC+1 

EMAX = 0. DO 

DC 997 I - 1 , NS 

VECP ( I , LC )-VECt I , K ) 

VECI< I , LC ) = V£C( I * L K ) 

CNOD=VECR M ,LC)**2+VECI U,LC>**2 
IF (CM0D-EMAX)997, 990,990 

990 EMAX = CMOD 
M= I 

99? CONTINUE 



18 53* 


VMR = VECR(M,LC) 

1854, 


VMI = V EC I ( M »LC 1 

1855. 


SMXIN=1.D0/E MAX 

18 56. 


DC 980 1=1. NS 

1 857. 


VR = VECR (I » LC) 

1858. 


VI = VECI (I » LC > 

1850. 


VECRN ( I .LC! =EMX IN4 (VR4VMR+V I*VM I ) 

I860. 

580 

VEC INI I f LC >=EMX IN* (-VR*VMI+VI *VMR | 

1861. 


VRVtLC l=WZC K) 

1862. 


VIV(LC>=WY(K> 

1863. 

592 

L M= 0 

1864. 


IF(0ABS(WY< K)+WY(LK) ).LT.1.C-10) IM1 

1865. 


GO TO 999 

18 66. 

558 

LR=LR+1 

1867. 


00 996 1=1, NS 

1868. 


RVECi I,LR)=VEC( I,K) 

1869. 

596 

REMOQ=RVEC( I *LR )**2+REMC0 

1870. 


RMOD=CSCRT( REMOD) 

1871. 


DC 995 1=1, NS 

1872. 

995 

R VEC ( t , LR )=R VEC ( I » LR ) /R MOD 

18 73. 


VRR V ( t R ) =W£ ( K ) 

1874. 


VRIV(LR)=WY(K) 

1875. 


IF(K.EC.NS) GO TO 1000 

1876. 

599 

CONTINUE 

1877. 

C1000 

WRITE C6.7) 

1878* 

1000 

CONTINUE 

1379. 


IF (LC .EQ. 0) GO TO 961 

1880. 


DC 960 J=l, LC 

1881. 


WRITE (6 ,101 J,J ,VRV( J 1, VI V( J» ,VECRN(1,J 1 ,VEC!NU, J) 

1882. 

960 

WR I TE t 6 ,1 1 1 (VECRMK.J) ,VECIMtK,J l,K=2,NS ) 

1833. 

961 

IF (LR.EQ.O) GO TO 1005 

1884. 


WR I TE ( 6 . 7 ) 

1885. 


DO 965 J= 1 , LR 

1886. 


WRITE16.R) J,J,VRRV(J),VRIV(J) ,RVEC(1,J ) 

1 887. 

965 

wR I TE ( 6 ,9 > ( RVEC( K ,J ) ,K=2 ,NS> 

13 88* 

1 Cu5 

RETURN 
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1 890* 
1891 . 
1892. 

18 93. 
1694 , 

1895. 

1896. 

1897. 

1898. 

1899. 

1900. 
1901 . 

19 02 . 

1903. 

1904. 
19C5. 

1906. 

1907. 
1 908. 

1909. 

1910. 

1911. 

1912. 

1913. 

1914. 

1915. 

1916. 
19 17. 

1918. 

1919. 

1920. 

1921. 

1922. 

1923. 

1924. 


c NL) 

SUBPOUT 1 NE GMADDINM, A,B,P,N ,M > 

REAl*8 A,B,R 

DIMENSION A ( NM) , B ( NM ) »R ( N M > 

DC 10 1=1, NM 
10 R(I > = A ( H*B( I) 

P FT URN 
END 

SUBROUTINE G*PRDl NM,ML,NL,A,B ,fi,N,M,L) 

P FAL* 6 A f B » P 

DIMENSION A< NM) ,3< ML) ,R(NL) 

IR = 0 
IK=-M 

DC 10 K=1,L 
IKMK+* 

OC 10 J = 1 , N 

I R = I R + 1 

J 1 = J-N 
I 6= I K 
R ( I R ) =0 . 

DC 10 1=1, M 
J I = J I + N 
I B= I B+ 1 

10 RII R) =R(I Rl +AIJ I )*B< I B) 

RETURN 

END 

SUBROUTINE GAI N < M ,NS ,NC , RM ,WR , K I , *F , GN , W 1 1 ,TC B , 
1W21,AT,LT,C,CI,CT,MHS ,MTl 
IMPLICIT REALMS <A-H,C-Z) 

D I MENS I CN WRIMJ ,WI <M) ,VF(M,M),GMNS,NS) ,RM(M,M) 

DIMENSION Will NS, NS) ,TCB< M ,NS> ,W21 (NS, NS ) , LT ( MHS ) , MT ( MHS ) 
0 1 MENS 1 ON AT (MHS ) 

DIMENSION C ( NSI ♦ C HNS » , CT I N S , NS > 

MF = NS 
K = 1 
KP = 1 



1925. 

1926. 

1927. 

1928 . 

1929 . 
19 30 . 

1931 . 

1932. 

1933. 
19 34. 

1935. 

1936. 

1937. 

1938. 

1939. 

1940. 

1941. 
19 42. 
1943. 
19 44. 

1945. 

1946. 

1947. 

1948. 

1949 . 
19 50 . 

1951 . 

1952 . 

1953. 

1954. 

1955. 
1956* 

1957. 

1958. 

1959. 

1960. 


KM = 1 

10 IF(K.GT.M) GO TO 200 

WZ=-WR(K)**2-WI IK M*2 4l . 

I FI WZ > 100,5000,50 
50 I F( WI (K ) I 80,75,80 

C EIGENVECTOR FOP REAL EIGENVALUE, POSITIVE 

75 CONTINUE 
KP = K P+1 
K=K + 1 
GO TO 10 

C EIGENVECTOR FOR COMPLEX E I GENVAL UE , PCS IT I VE REAL PART 

80 CONTINUE 

KP = KP+2 , 

K = K + 2 
GO TO 10 

100 IF ( W I ( K > ) 120,110,120 

C EIGENVECTOR FOR -REAL E I GENVALUF , NEGATI VE REAL PART 

110 C ( KM ) = WR<K> 

CI(KN) = MKK) 

DO 95 J= 1 ,M 
95 TCB( J , KM) = VF{ J , K ) 

KN = KN+1 
K=K + 1 
GO TO 10 

C EIGENVECTOR FOR COMPLEX El GENVAL UE , NEGATI VE REAL PART 

120 RR = WRCK) 

R I = MIIKI 
CfKK) = RR 
CIKN+1) = RR 
CHKN) = Rl 

Cl C KN + 1 ) = - RI 
DO 121 J = 1 ,M 
FP = VF ( J , K I 
FI = VF{ J ,K + 1 ) 

TCB(J»KN) = FR+FI 

121 TCB(J,KN+1 > = FR-FI 
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19 61, 
1962. 
19 63 . 
1 96*. 
196?. 
1966. 
19 67 . 

1968. 

1969. 

1970. 
19 71. 
1972. 
1 9 73. 
19 7*. 
1979. 

1976. 

1977. 

1978. 

1979. 
19 80. 

1981. 

1982. 

1983. 
1 9 89. 
1985. 
19 86 . 
19 87. 

1988. 

1989. 

1990. 

1991. 

1992. 

1993. 

1994. 

1995. 
19 56, 


KN = KN + ? 

K - K + 2 
GO TO 10 
200 CONTINUE 

DO 29° I = 1 * NS 
00 299 J = 1 . NS 
295 CT( I* J> = TCB{ 1 ,J) 

C FORMATION CF Wll 

DO 300 I * 1*MH 
cn 300 J = 1,MH 
3 CO Wild ,J> = TC3( I * J) 

KNS = NS+1 
C FORMATION OF w21 
DO 320 1=1, M H 

00 32.0 J=1,MH 

320 w21 ( I ♦ J )= TC8( ! +MH »J) 

322 CONTINUE 
MHS=MH*MH 
C INVERT Wll 

DC 327 J=1,NS 
DC 327 K- 1 , NS 
I=( J-l l*NS +K 
327 AT( I)=W11(K, J) 

CALL MINVtMHS, AT , MH, DE TC , L T, MT ) 

DC 328 J-1,NS 
DC 328 K=1»NS 

1 = < J— 1 > *N S + K 
3 28 WlltK* J)=AT(I) 

C 

C CALCULATE THE RGAIN MATRIX 
DO 325 IL=1,MH 
DO 325 JL=1,MH 
GMILfJL) * O.DO 
DO 325 KL=l f MH 

32 5 GNC IL,JL) = GNIIL ,JL 1 + W 21 ( I L , KL ) *Wl 1 ( KL » J L ) 

GO TO 6000 
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i. ''UOU WM Itlrdl 

l'JSe. 1 ^CRWAT( «1 • ,• 7ERC EIGENVALUE -RICATT I -QUATICN hA$ NH SSM 

199=5* 6000 RETURN 

2000. END 

20C1. /* 

2C02. //LKED. SYSLMOC OD ESN = $ ?70. ER AM 0 ISC ) , 

2003 . // UNIT=231<t t VDL = SER=SY$14 f CI$P=< t KEEFI , 

2-004- • // SPACE-ITRK »( 5, 5,i> ,RLSE) 

20C5. /* 

l 




Appendix B 


THE SIMULATION SCHEME : A SIMULATION ALGORITHM FOR A DIS - 

CRETELY CONTROLLED, LINEAR, CONTINUOUS, TIME INVARIANT SYSTEM 


If one is interested in investigating the behavior between the sampl 

ings, say, n times during the interval, T, the basic integration 
S 

step will be T (T = T/n ). The integration algorithm solves, step 
s s s 

by step, the following equations. 


x J , = ® (T )x. + I\(T )u.. 

i+1 P x s' 1 P s' j 


(B-l) 


The control u. is constant during n iterations. $ and T are 
J S F Ir 

calculated a priori from 


<t> (T ) = e 

P^ s' 


r r,( T ) = 

S 


" F P T s 


f s .'V 


dr 


(B-2) 


% * 


u. is calculated every T (after n iterations of Eq, B-l) from: 
3 s 


u . = Cx . + u(T) 

3 3 ° 

x , = (I - KH)[($(T) + r(T)*C)x. + u(t)J + Ky. 


j+1 


i+1 


(B-3) 


Vl * "V 


Equations B-3 are the compensator equations. <D and P are calculated 
similarly to (B-2). Equations B-3 are mechanized in the discrete com- 
pensator, the derivation of the equations having been given in Chapter 
VI. 
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All the matrices of (B-3) are calculated by our discrete synthesis 
program, DISC. $ and F represent the assumed system, while and 
belong to the actual real plant. The distinction between these two 
systems is necessary in order to investigate the sensitivity (Ch. Vi). 

The relation between T, T , and n is described in Fig. B-l. 

s s 



FIG, B-l THE SIMULATION TIMING* 

T - sampling time of the plant. 

T = sampling step in the simu— 
s i . i 

lat ion. 
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